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INTRODUCTION 


Two  phase  gas-liquid  flow  and  its  associated  interfaces  exist  in 
a  wide  variety  of  situations  of  importance  to  the  Navy  and  this  has 
prompted  the  study  of  the  basic  flow  mechanics  which  underlie  this 
complex  process.  The  existence  of  wind-wave  interactions  over  large 
bodies  of  water  have  long  been  recognized  as  a  special  case  of  two 
phase  flow  where  the  presence  of  the  deformable  interface  plays  a 
commplex  role  in  the  generation  of  waves  due  to  the  action  of  the 
wind.  Less  well  recognized,  but  of  great  importance,  are  situations  of 
two  phase  flow  which  are  found  in  componenet  of  power  systems 
such  as  condensedrs,  boilers  refrigeration  loops  and  cryogen  lines. 
Here  the  characteristics  of  two  phase  flow  are  critical  to  the  reliable 
design  and  safe  operation  of  such  systems.  A  basic  understanding  of 
gas-liquid  flows  is  also  central  to  the  design  of  safety  systems  for 
nuclear  reactors.  In  fact  some  of  the  earliest  research  contributions 
to  the  mechanics  of  two  phase  flow  was  motivated  by  the  need  to 
design  such  emergency  safety  systems  for  the  early  versions  of  the 
Navy  nuclear  submarine  fleet. 

During  simultaneous  flow  of  gas  and  liquid  in  a  conduit  the 
phases  distribute  in  a  variety  of  patterns.  These  include  gas  bubbleds 
distributed  in  the  liquid,  alternating  flows  of  gas  and  liquid  plugs  and 
annular  flow  where  the  two  phases  are  separated  with  the  liquid 
flowing  predominately  along  the  wall  and  the  gas  flowing  in  the  core. 
It  is  the  latter  configuration  which  has  been  the  subject  of  this  study. 
In  this  configuration  the  interface  is  covered  by  a  complex  chaotic 
system  of  waves  of  considerable  ampolitude  compared  to  the  thin 
film  over  which  it  rides.  Furthermore  the  high  speed  gas  causes 
detachment  of  the  some  of  the  liquid  from  the  waves  and  this  forms 
an  entrained  droplet  phase  which  travels  with  the  gas.  IN  this 
research  program  this  complec  problem  has  been  attacked  from  a 
number  of  directions. 

A.  Phase  plane  and  bifurcation  analysis  of  the  complex  wavy 
structure  of  the  interface.  These  results  appeared  in  two  published 
papers  attached: 

•  Phase  plane  and  bifurcation  analysis  of  thin  wavy  films 

under  shear.  AIChE  J  15^177-186  (1989) 


•  Methods  of  deterministic  chaos  applied  to  flow  of  thin 
wavy  films.  AlChE  J  22  481-489  (1991) 

B.  Experimental  and  numerical  investigations  of  the  mechanics 
of  wavy  flow  including  the  solution  of  the  momentum  equation  for 
the  velocity  field  under  the  free  wavy  surface  and  the  interactions 
between  the  waves  and  this  field.  These  appear  in  two  publications 
which  are  attached. 

•Insights  into  the  hydrodynamics  of  free  falling  wavy 
film.  AICHE  J.  22.187-195  (1989) 

•  Numerical  investigation  of  large  wave  interactions  on 
free  falling  films.  Int.  J.  Mult.  Flows  12,  357-370  (1989) 

C.  Studies  of  mass  transfer  and  the  role  of  the  waves  on  the 
enhancement  of  the  transfer  process.  These  work  is  presented  in 
two  publications  attached. 

•  A  numerical  study  of  mass  transfer  in  free  falling  wavy 
film.  AIChE  J  1379-1390  (1990) 

•An  experimental  study  of  mass  transfer  from  a  wall  into 
a  wavy  falling  film.  Chem  Eng  Sci  42  4323-4331  (1992) 

Two  studies  initiated  and  carried  out  under  this  study  and  now 
nearing  completion  are 

•  Studies  of  turbulent  gas  flow  over  a  wavy  interface. 

•  Entrainment  from  a  wavy  interface  and  drop  dynamics 
of  the  entrained  phase. 

PhD  theses  for  each  of  these  areas  are  now  essentially  complete. 
Papers  have  been  accepted  for  presentation  on  both  subjects  at  the 
next  annual  meeting  of  the  American  Institute  of  Chemical  Engineers. 
Copies  of  the  manuscripts  will  be  forwarded  when  ready.  It  is 
estimated  that  they  will  be  available  about  the  end  of  the  year. 


Methods  of  Deterministic  Chaos  Applied  to  the 
Flow  of  Thin  Wavy  Films 

C.  1'..  Lacy.  M.  ShciiUuch,  and  A.  E.  Pukler 
Dept,  of  Chemical  Engineering,  University  of  Houston,  Houston,  IN  77204 

The  structure  of  thin,  wavy  falling  films  was  studied  to  evaluate  whether  the 
random-appearing  wave  structure  is  a  result  of  deterministic  chaos  or  a  purely 
stochastic  process.  The  lime-varying  film  thickness  was  obtained  at  different  spatial 
locations  near  the  point  of  wave  inception  for  flow  rates  in  the  range  of  Re  = 3-10 . 

Under  all  conditions  the  wave  structure  was  aperiodic  in  nature  and  displayed  none 
of  the  known  transitions  to  chaos.  However,  the  power  spectra  followed  an  expo¬ 
nential  decay  law  at  high  frequencies  that  is  characteristic  of  chaotic  systems.  The 
estimated  attractor  dimension,  used  to  characterize  the  complexity  of  a  chaotic 
system,  was  much  higher  than  those  of  known  model  chaotic  systems.  It  is  dem¬ 
onstrated  that  these  high  values  could  be  explained  due  to  small  levels  of  noise 
present  in  experimental  situations.  Since  experimental  data  are  seldom  noise  free, 
a  basic  limitation  in  applying  these  methods  to  experimental  measurements  is  dem¬ 
onstrated. 


Introduction 

The  hydrodynamic  behavior  of  thin  wavy  falling  films  has 
been  a  subject  of  intensive  investigation  for  about  forty  years. 
These  films  are  widely  employed  in  equipment  for  heat  trans¬ 
fer.  mass  transfer,  and  chemical  reacting  systems.  In  addition 
to  the  practical  need  for  understanding  the  mechanics  of  this 
type  of  flow,  there  are  challenging  theoretical  problems  em¬ 
bedded  in  the  task  of  modeling  these  wavy  films.  This  com¬ 
bination  has  given  rise  to  an  extensive  literature  on  this  subject. 

The  earliest  work  was  based  on  the  use  of  integral  equations 
of  the  boundary  layer  type  to  solve  the  equations  of  motion 
(Kapitza  and  Kapitza,  1949;  Shkadov,  1967).  These  approaches 
were  based  on  the  assumed  existence  of  a  periodic  interface 
and  produced  first-order  estimates  for  wavelength  and  velocity. 
A  long  series  of  papers  including  work  by  Benney  (1966),  Lin 
(1969),  and  Whitaker  (1964),  used  linear  stability  analysis  to 
find  the  wavelength  and  velocity  of  the  fastest  growing  wave, 
again  assuming  a  periodic  perturbation.  This  initial  periodic 
disturbance  is  thought  to  evolve  into  the  more  complex 
waveforms  observed  in  experiments  as  a  result  of  the  nonlinear 
nature  of  the  equations.  Nonlinear  stability  analyses  have  also 
been  pursued  (Pumir,  1983).  Recent  work  on  the  nonlinear 
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nature  of  wavy  films  (Chang,  1987;  Sheintuch  and  Dukler, 
1989)  searched  for  infinitely  long  periodic  waves  by  finding 
the  conditions  for  existence  of  the  homoclinic  orbit. 

But  even  a  cursory  examination  of  measured  wave  traces 
raises  some  doubts  as  to  the  usefulness  of  the  idea  of  a  small- 
amplitude  periodic  wave  as  the  model  for  the  initial  phase  of 
the  wave  motion  or  of  isolated  waves  as  a  model  for  the  de¬ 
veloped  ones.  Figure  Id  shows  a  wave  trace  for  a  falling  liquid 
film  of  water-glycerine  solution  taken  with  a  conductivity  probe 
mounted  in  a  vertical  pipe  of  50.8  mm  dia.  as  described  below. 
The  flow  rate  corresponds  to  a  Reynolds  number  of  3.9.  The 
probe  was  located  0.346  m  below  a  carefully  leveled,  sharp- 
edged  overflow  weir  that  served  as  the  feed  device.  The  film 
thickness  data  are  shown  after  low-pass  digital  filtering  at  25 
Hz  to  remove  noise.  At  this  location  the  wave  amplitude  is 
less  than  0.2S^b  of  the  mean  film  thickness.  At  all  positions 
closer  to  the  feed  the  waves  were  so  small  that  they  could 
barely  be  detected  even  with  the  special  circuitry  used  for  this 
purpose.  Note  that  while  the  period  between  successive  waves 
is  quite  regular,  the  amplitude  is  very  tandom.  Kapitza  and 
Kapitza  (1949)  in  their  classical  study  of  waves  on  falling  films 
found  it  necessary  to  pulse  the  feed  to  produce  periodic  waves. 
In  the  absence  of  pulsing  they  too  reported  that  the  waves  were 
random  in  appearance.  Thus  one  must  question  whether  the 
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c  onsider  a  set  of  <V  points,  constructed  as  in  Liq  I.  that  lie 
on  the  reconstructed  attractor.  Arbitrarily  select  a  point  x  tin 
this  attractor  as  a  reference  point.  Now  choose  at  random  a 

subset  of  A'  points  denoted  by  y,C  =  1.2 . A-  and  k<N)  from 

the  original  set  of  A'  points  and  consider  the  distance  from  x 
to  each  point  y,.  Define  6  as  the  single  minimum  of  these 
distances,  that  is,  the  distance  to  the  nearest  neighbor. 

6  =  minlx-yfl  (2) 

To  obtain  a  statistically  useful  value  of  the  minimum  distance, 
this  calculation  is  repeated  over  many  randomly  chosen  ref¬ 
erence  points  and  an  average  obtained.  (6).  The  process  is 
then  repeated  for  a  sequence  of  k  values  up  to  k  -  N~  I  for 
each  x.  The  number  of  nearest  neighbors  contained  in  a  D- 
dimensionai  hypersphere  of  radius  5  around  a  given  point 
should  vary  as  if  the  attractor  is  cf-dimensional.  Thus  it  is 
argued  that 

<6> (3) 

Hence 

log<6> — — -log  k  (4) 

a 


Figure  3.  Low  Re  flow  system. 


ing  values  of  r.  The  optimal  reconstruction  occurs  for  the 
smallest  value  of  r  for  which  x(r)  and  x(r-t-T)  can  first  be 
considered  independent.  The  criterion  for  choosing  this  value 
of  r  is  the  first  minimum  in  the  index  of  mutual  information, 
/,  as  described  by  Fraser  and  Swinney  (1986).  This  index  is  a 
measure  of  the  degree  of  predictability  of  a  measurement 
x (t  +  r)  given  a  measurement  x (/). 


The  negative,  inverse  slope  of  a  log  <6)  vs.  log  k  plot  is  the 
fractal  dimension. 

All  experimental  measurements  of  forced  flow  systems  in- 
lude  electronic  noise  as  well  as  noise  due  to  random  vibrations 
reaching  the  system  through  piping  or  pumps.  Small-amplitude 
noise  can  be  expected  to  distort  distances  between  closest  near- 
si  neighbors  and  their  reference  point.  To  partly  alleviate  this 
rror  Kostelich  and  Swinney  (1987)  suggest  that  6  be  calculated 
for  the  tenth  or  hundredth  nearest  neighbor.  For  each  choice 
'  D  and  nearest  neighbor,  it,  an  estimate  of  d  can  be  calculated 
turning  k  is  kept  sufficiently  larger  than  it.  The  value  of  d 
is  said  to  converge  to  the  attractor  dimension  when  d  becomes 
independent  of  D  and  n  as  both  are  increased.  Additional 
pects  of  the  noise  problem  are  discussed  later  in  this  paper. 
The  computed  value  of  d  has  been  shown  to  depend  on  the 
value  of  the  delay  time,  t,  used  to  construct  the  D-dimensional 
vectors  in  Eq.  1.  If  r  is  too  small,  then  each  x(r)  approaches 
/  +  r),  and  the  reconstructed  attractor  will  be  a  45*  plane  in 
:  ;  phase  space.  For  large  values  of  r,  the  attractor  dimension 
lends  to  approach  the  embedding  dimension,  D,  of  the  rccon- 
r  ucicd  phase  space  due  to  the  stretching  and  folding  nature 
t  chaotic  systems  (Fraser  and  Swinney,  1986). 

fo  obtain  the  optimal  delay,  a  series  of  two-dimensional 
reconstructions  of  |x(r).  *(/  +  r)  |  arc  generated  with  incrcas- 


Experimental  Equipment 

Figure  3  is  a  diagram  of  the  experimental  flow  loop  used 
for  studies  of  free-falling  liquid  films  at  low  Re{Re-  4Q/v 
=  3-10).  Where  Q  is  the  volumetric  flow  rate  per  unit  perimeter 
and  v  is  the  kinematic  viscosity.  The  measurement  section 
consisted  of  a  50.8  mm  dia.  tube  0.47  m  long.  The  liquid  feed 
tank  contained  a  sharp-edged  circular  weir  over  which  the 
liquid  flowed  from  the  reservoir  into  the  pipe.  The  distance 
from  the  sharp  edge  to  the  bottom  of  the  flange  was  86  mm. 
This  was  followed  by  a  development  section  having  one  of 
three  lengths:  86, 138,  or  284  mm.  Wave  motion  is  not  observed 
immediately  after  the  film  is  formed  at  the  overflow.  The 
development  section  provided  the  length  needed  at  different 
flow  rates  to  cause  the  waves  to  first  appear  in  the  measuring 
section  where  film  thickness  probes  were  located 

I nstanianeous  film  thickness  data  were  obtained  using  closely 
spaced  parallel  wires  to  make  resistivity  measurements,  which 
are  related  to  film  thickness  by  calibration.  Details  of  this 
method  including  the  electronic  circuit  used  appear  in  the  thesis 
by  Zabaras  (1985).  For  these  low  Re  measurements,  an  array 
of  seven  probes  was  spaced  in  the  direction  of  flow  with  dis¬ 
tances  from  the  feed  point  listed  in  Tabic  I .  A  glycerinc-watcr  - 
NaOH  solution  having  a  kinematic  viscosity  of  5.6x10 
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Figure  5.  Experimental  film  thickness  traces  at  Re  -  7.5. 

Distance  belo*  overflow  weir: 

a.  0.21  m;  b.  0.26  m,  c.  0.41  m;  d.  0.61  m 

of  all  the  waves  returning  to  a  single  thickness  and  the  variation 
of  wave  amplitude  being  reflected  in  differences  in  the  peak 
value.  Figure  7c  shows  this  structure  in  which  the  lefihand 
comer  resembles  motion  in  the  vicinity  of  a  saddle  point.  With 
further  distance  in  the  axial  direction  the  wave  periods,  which 
heretofore  had  been  regular,  now  become  chaotic,  as  shown 
in  the  time  trace,  the  spectrum,  and  in  Figure  7d. 

When  two  spectral  peaks  of  about  equal  power  exist  at  wave 
inception  (Re=  5.6),  the  approach  to  a  more  uniform  wave 
amplitude  is  not  observed  and  the  waves  progress  more  rapidly 
to  the  condition  of  chaotic  wave  periods. 

The  power  spectra  provide  a  suggestion  that  the  underlying 
process  may  be  a  manifestation  of  deterministic  chaos.  Sigeti 
and  Horsthemke  (1987)  have  shown  that  for  a  yth-order  sto- 
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Figure  6.  Power  spectral  density  of  Re  =  7.5  data  in  Fig¬ 
ure  5. 


Film  Thickness  (mm) 

Figure  7.  Phase  plane  portraits  of  Re  =  7.5  data  in  Figure 
5. 

chasiic  differential  equation  the  spectral  power  at  high  fre¬ 
quencies  follows  a  power  law  decay,  S  =  cf~:'.  Deterministic 
equations  are,  on  the  other  hand,  infinitely  differentiable  and 
the  spectrum  must  decay  exponentially  at  hign  frequencies.  All 
spectra  observed  were  characterized  by  exponential  decay. 

Attractor  Dimension.  For  each  data  set  listed  in  Table  1, 
the  nearest-neighbor  algorithm  for  the  computation  of  the 
attractor  dimension  was  applied  using  the  following  parame¬ 
ters: 


Total  points,  N  =  100,000 

Embedding  dimension.  D  =  8-20 

Nearest  neighbors,  n  =  100-600 

Subsets  of  N  points,  k  =  15,000-100,000 


These  parameters  resulted  in  processing  2.000-4,000  waves 
with  each  wave  represented  by  30-40  discrete  points. 

Figure  8  illustrates  the  method  of  Badii  and  Poiili  (1985) 
when  applied  to  the  time  series  data  from  probe  5  at  Re  =  7.5. 
Figure  8a  is  a  plot  of  log2  <6>  vs.  log;  k  for  one  D  and  n  pair. 
An  estimate  of  the  attractor  dimension,  </,  is  obtained  from 
the  slope  as  given  by  Eq.  4.  The  results  of  repeating  this  process 
for  each  combination  of  D  and  n  appear  in  Figure  8b.  Here 
the  dimension  has  converged  with  increasing  D  and  n.  indi¬ 
cating  an  attractor  dimension  of  about  3.5.  There  are  instances 
where  convergence  was  not  observed  with  increasing  D ,  as 
illustrated  in  Figure  9,  although  convergence  with  increasing 
n  was  always  observed.  However,  this  system  docs  not  rep¬ 
resent  purely  stochastic  noise,  since  under  that  condition  one 
would  expect  d  to  be  close  to  D.  This  behavior  seems  to  be 
related  to  the  presence  of  noise  in  the  signal  and  presents  a 
basic  problem  in  processing  and  analyzing  data  from  real  sys¬ 
tems  as  compared  to  analyzing  data  from  mathematical  models 

Figure  10  shows  the  trends  in  complexity  of  the  time  traces 
as  position  and  flow  rates  are  changed.  Estimated  dimensions 
appear  here  that  have  been  computed  for  D  -  20,  N  =  100.000. 
n  =  600  using  1 .000  reference  points  In  some  cases  definitive 


AlCtiK  Journal 


April  I'W I  Vol.  37,  No.  4 


485 


Amplitude 


i 


J _ I - 1 - 1 _ 1 _ 1 _ L 


Value  of  j 

Figure  11.  Time  trace  fora  sine  wave  with  multiplicative 
noise. 

a.  *  =  0.125.  o.  *  =  0.25 
c.  *  =  0.5;  d.  *  =  1.0 


Attempts  to  diminish  the  effects  of  noise  by  low-pass  filtering 
had  negligible  impact  on  the  calculated  dimension  of  the  at¬ 
tractor.  Methods  for  eliminating  noise  have  recently  been  ad¬ 
vanced  by  Kostelich  and  Yorke  (1988),  but  they  appear 
applicable  only  when  the  noise  is  small  (<!%),  a  condition 
that  does  not  appear  to  exist  here.  Until  improved  methods 
are  available,  it  would  appear  that  this  method  for  determining 


Value  of  j 

Figure  1 2.  Time  trace  for  a  sine  wave  with  additive  noise. 

a  *  =  0  025:  b  *  =  0  15.  c.  *  =  0  5 


Embedding  Dimension,  D 
Figure  13.  Dimension  estimates  for  a  sine  wave  with 
multiplicative  noise. 

*.*  =  0.125;  b.  *=0.25 
c.  *«0.5;  d.  *=  1.0 


if  a  process  displays  the  characteristics  of  deterministic  chaos 
will  be  of  very  limited  utility  when  applied  to  experimental 
data  from  physical  systems. 


High  Reynolds  number  flows 

Figure  2  is  representative  of  the  time  trace  obtained  at  higher 
flow  rates.  The  spectrum  appears  in  Figure  15.  Attractor  di¬ 
mensions  were  calculated  for  ten  conditions  in  the  range 
Re=  310-3,100  and  include  data  for  free-falling  films  as  well 
as  those  with  countercurrent  and  cocurrent  intcrfacial  shear. 
All  of  these  traces  were  characterized  by  poor  convergence  and 
high  dimension.  Extensive  studies  did  not  reveal  any  coherent 
dependence  on  Re  or  on  the  amount  or  direction  of  the  in¬ 
terfacial  shear.  However,  the  high-frequency  end  of  the  spec¬ 
trum  still  indicates  an  exponential  decay,  suggesting 
deterministic  chaos.  It  is  likely  that  noise  plays  a  large  role  at 
these  high  rates  where  isolation  of  the  system  is  more  ifficult. 


Discussion 

An  extensive  literature  exists  for  analyzing  time-dependent 
data  using  methods  of  fractal  geometry.  When  the  data  arc 
developed  from  mathematical  models  these  methods  of  anal¬ 
ysis  provide  new  insights  into  the  behavior  of  these  nonlinear 
equations.  However,  attempts  to  analyze  data  obtained  dircctlv 
from  experiment,  as  has  been  presented  in  this  paper,  face 
severe  difficulties.  Such  data  arc  accompanied  by  significant 
amounts  of  noise,  and  it  is  shown  that  noise  levels  common!- 
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A  Numerical  Study  of  Mass  Transfer  in  Free 
Falling  Wavy  Films 


Numerical  simulations  of  mass  transfer  into  falling  liquid  films,  both 
through  the  wavy  interface  and  from  the  wall,  have  been  performed  for 
experimentally  measured  large  waves  within  which  the  flow  fields  have 
been  computed.  Experiments  have  shown  that  the  occurrence  of  waves 
on  free  falling  films  causes  dramatic  increases  in  mass  transfer  into  the 
film,  even  under  laminar  flow  conditions.  Wave  effects  have  been 
modeled  in  several  ways,  none  of  which  predicts  the  observed  rate  of 
enhancement.  The  present  numerical  procedure  includes  solving  the 
convective-diffusion  equation  for  wavy  films  by  extending  a  technique 
developed  for  hydrodynamic  simulation.  The  presence  of  waves  is 
shown  to  cause  significant  velocities  normal  to  each  interface.  In 
conjunction  with  recirculation  within  the  large  waves,  these  flow  pat¬ 
terns  produce  transfer  rates  for  large  waves  that  are  several  times 
larger  than  predicted  for  quasiparallel  velocity  fields.  Experimental 
wave  structure  data  were  used  to  define  the  dimensions  and  frequency 
of  an  average  large  wave  and  surrounding  substrate.  Computed  transfer 
rates  at  both  the  gas-liquid  interface  and  the  wall  for  a  film  composed  of 
a  periodic  sequence  of  average  waves  agree  well  with  published  data. 
These  simulations  confirm  the  inadequacy  of  parabolic,  or  Kapitza-type 
velocity  profiles  in  formulating  transport  models. 


Frederic  K.  Wasden 
A.  E~  Dukler 

Department  ot  Chemical  Engineering 
University  of  Houston 
Houston.  TX  77204 


Introduction 

The  ability  of  liquid  films  to  transfer  large  amounts  of  heat  or 
mass  with  low  hydraulic  resistance  has  led  to  their  use  in  a  wide 
variety  of  industrial  processes.  For  unit  operations  in  which  the 
liquid  phase  is  distributed  as  a  film,  the  throughput  of  the  unit  is 
often  determined  by  the  liquid-phase  transport  resistance; 
trickle-bed  and  failing-film  reactors,  wetted-wall  absorbers,  and 
vertical  evaporators  are  examples.  Upon  contacting  a  solid 
surface,  liquid  films  quickly  evolve  to  a  complex  array  of  waves 
whose  amplitudes  vary  from  much  less  to  much  greater  than  the 
mean  thickness.  This  gravity-driven  behavior  is  observed  for  all 
flow  rates  of  industrial  importance,  even  in  the  absence  of 
interfacial  stresses  due  to  adjacent  gas  flow.  Figure  1  shows  a 
sample  lime  tracing  of  the  interface  of  a  fully  developed  laminar 
film  falling  freely  down  a  vertical  tube  without  gas  flow.  The 
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presence  of  large  waves  flowing  over .'  thin  substrate  which  itself 
is  covered  by  small  waves  is  clearly  displayed. 

Mass  transfer  in  liquids  is  characterized  by  extremely  long 
time  scales  for  molecular  diffusion.  The  ratio  of  thermal  to 
molecular  diffusivities  in  liquids  is  generally  greater  than  100, 
suggesting  that  mass  transfer  rates  reflect  fluctuations  in  flow 
field  to  a  greater  extent  than  do  heat  transfer  rates.  This 
speculation  has  been  confirmed  by  reviews  of  experimental 
studies  (Seban  and  Faghri,  1978;  Henstock  and  Hanratty, 
1979).  This  disparity  is  exploited  in  the  present  work  to  provide 
a  comprehensive  examination  of  a  combination  of  numerical 
procedures. 

A  simplified  view  of  a  liquid  film  suggests  that  transport  is 
limited  by  diffusion  in  the  direction  normal  to  the  transfer 
surface.  The  velocity  of  a  contaminant  traveling  normal  to  the 
interface  as  the  result  of  diffusion  is  of  the  order  of  the  ratio  of 
the  diffusivity  to  the  film  thickness.  The  comparable  velocity  in 
the  streamwise  direction  due  to  advection  is  of  the  order  of  the 
average  film  velocity.  The  ratio  of  the  advective  to  the  diffusive 
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f  .cjjre  1.  Film  thickness  time  trace,  Re  —  880. 


velocity,  the  Peclet  number,  is  of  order  1  O'*—  1 0*  for  industrially 
important  thin  film  flows.  Experiments  have  shown  that  the 
presence  of  waves  on  films  causes  dramatic  increases  in  heat  or 
mass  transfer,  even  for  laminar  flows.  While  a  wide  range  of 
wave  amplitudes  exists,  it  is  speculated  tha‘  the  large  waves, 
ranging  from  two  to  five  times  the  substrate  thickness  in 
amplitude,  control  the  transport  rates  (Dukler,  1977). 

Predicting  the  enhancement  of  transport  rates  due  to  the  wavy 
interface  has  provoked  many  studies  of  film  hydrodynamics. 
Since  first  being  addressed  by  Kapitza  and  Kapitza  (1949), 
studies  of  the  linearized  hydrodynamics  within  wavy  films  have 
yet  to  predict  the  wide  variety  of  wave  shapes,  sizes,  and  speeds 
observed  experimentally.  Numerical  studies  of  the  problem  are 
limited.  Bach  and  Villadsen  (1984)  succeeded  in  predicting 
velocity  profiles  in  traveling  waves  using  a  finite-element  tech¬ 
nique.  but  their  results  were  limited  to  Reynolds  numbers  less 
than  100.  The  film  Reynolds  number  is  defined  as  Re  =  4 Q/v, 
where  Q  is  the  mass  flow  rate  per  unit  perimeter  and  v  is  the 
kinematic  fluid  viscosity.  Recent  numerical  studies  of  hydrody¬ 
namics  in  isolated  and  interacting  large  waves  at  a  Reynolds 
number  of  880  (Wasdcn  and  Dukler,  1989a,  b)  have  provided 
information  on  the  flow  fields  existing  in  these  waves.  The 
hydrodynamic  studif  predict  regions  of  large  streamwise  accel¬ 
eration  in  the  waves  as  well  as  recirculating  zones.  Nakaya 
( 1 989)  has  also  found  these  flow  patterns  in  large  waves  at  low 
flow  rates. 

This  study  focuses  on  numerical  prediction  of  mass  transfer  to 
laminar  free  falling  films  for  both  gas  absorption  at  the  free 
surface  and  dissolution  at  the  wall  in  the  presence  of  compli¬ 
cated  velocity  fields.  The  former  situation  has  been  studied 
extensively,  both  experimentally  and  theoretically,  and  is  re¬ 
viewed  by  Henstocx  and  HanraUy  (1979).  Information  u.n  the 
wall  transfer  is  much  more  limited. 


Figure  2.  Absorption  through  the  surface  of .  flat  film. 
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where  c(x,  y,  t)  is  the  concentration  of  the  species  and  u(x,  y,  t) 
and  »(x,  y,  r)  are  the  velocity  components  in  the  streamwise  (x) 
and  normal  (y)  directions.  For  steady  flat  film  flow,  the 
parabolic  velocity  profile  is  used  for  the  axial  velocity  and  the 
convective-diffusion  equation  reduces  to 
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where  streamwise  diffusion  is  negligible  at  the  high  Peclet 
numbers  experienced  in  practice.  For  absorption,  the  boundary 
and  initial  conditions  assume  an  initially  pure  solution  with  a 
saturated  interfacial  condition  and  no  flux  at  the  wall: 

c(x, y)  =  0  at  x  =  0  for  ally  (3a) 

c(x,  /i)  =  c„,uriltd  —  c,  at  y  =  A  for  all*  (3b) 

dc 

—  =  0  at  y  =  0  for  all  x  (3c) 


Existing  Models  and  Experimental  Studies 
Mass  transport  in  fiat  films 

Early  attempts  at  describing  transport  in  films  were  limited  to 
situations  where  the  film  was  assumed  to  be  flat.  Transfer  from 
the  adjacent  gas  phase  into  the  film  through  the  free  interface  is 
illustrated  in  Figure  2.  The  transport  is  governed  by  the 
convective-diffusion  equation,  which  can  be  written  in  two 


The  first  analytical  solution  of  Eqs.  2  and  3  was  presented  by 
Pigford  (1941),  who  assumed  that  the  gas-liquid  contact  time 
was  short  cnougii  that  variations  in  concentration  occur  only  in  a 
boundary  layer  near  the  free  interface.  Replacing  the  condition 
of  Eq.  3c  with  c(x,  0)  =  0  enabled  him  to  find  a  closed-form 
solution  for  the  local  mass  flux  at  any  position  x,  and  for  the 
mass  flux  averaged  over  a  column  length  of  L.  The  latter  solution 
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is  given  bclo*  and  is  referred  lo  as  the  short  contact  tune  theory 
(SCTT)  mass  transfer  rate, 


WSCT 


7DghH. 


(4) 


Mass  transfer  from  the  solid  wall  into  the  falling  liquid  film  is 
seldom  encountered  in  industrial  practice,  although  its  analog, 
the  process  of  heat  transfer  from  the  wall  is  quite  common. 
Figure  3  illustrates  the  mass  transfer  situation  for  a  fiat  film  in 
which  a  portion  of  the  wall  is  considered  to  be  active,  or  capable 
of  supplying  mass  to  the  film.  Boundary  and  initial  conditions 
are  similar  to  those  for  the  absorption  problem.  An  initially  pure 
solvent  is  assumed  to  be  saturated  at  the  liquid -solid  interface 
and  no  transport  is  allowed  at  the  free  interface: 


c(x,.y)  = 

0  at  x  =  0 

(3a) 

=  c,  for 

y  *=  0,  0  <  x  <  La 

(5a) 

3*1* 

II 

o 

at  y  *»  h 

(5b) 

0  at  y 

=  0  for  x  >  La 

(5c) 

No  general  analytical  solution  to  the  governing  equations 
exists  for  arbitrary  lengths  of  the  active  region,  LA.  For  small 
diffusivities  (large  Schmidt  numbers)  and  small  contact  times,  a 
boundary  layer  analysis  similar  to  Pigford’s  was  proposed  by 
Acrivos  (1960).  In  the  limit  of  low  Peclet  numbers  (high  mass 
dififusivity),  Spence  and  Brown  (1968)  solved  the  transport 
problem  using  a  Frobenius  scries  to  solve  the  ordinary  differen- 


Flflure  3.  Mas*  transfer  from  solid  boundary  into  flat  film. 


lial  equation  generated  by  the  Laplace  transform  of  the  original 
partial  differential  equation.  In  principle,  this  method  can  be 
extended  to  the  general  case,  but  ihe  Frobenius  series  docs  not 
possess  sufficiently  robust  convergence  characteristics  to  allow 
its  use. 

Gas  absorption  through  the  wavy  interface 

Most  gas  absorption  experiments  that  have  been  reported 
provide  information  on  mass  flux  or  mass  transfer  coefficients 
averaged  over  the  entire  length  of  film.  Simultaneous  local 
measurements  of  time-varying  wavy  film  thickness  and  concen¬ 
trations  have  not  been  reported,  so  little  information  is  available 
on  the  local  transport  process.  Emmert  and  Pigford  (1954) 
reported  mass  transfer  rates  in  agreement  with  short  contact 
time  predictions,  Eq.  4,  when  a  surfactant  was  added  to  suppress 
waves.  In  the  absence  of  surfactant,  mass  transfer  rates  two  (o 
three  times  greater  than  observed  for  fiat  films  were  measured 
for  Reynolds  numbers  of  200-1,200  and  Schmidt  numbers 
between  400  and  500  in  a  column  1.1m  long. 

Kafesjian  et  al.  (1961)  examined  the  rates  at  which  a  species 
is  absorbed  and  desorbed  from  a  film.  No  equivalent  to  the  short 
contact  time  theory  exists  for  desorption  since  the  concentration 
is  uniform  in  the  incoming  fluid.  Measured  values  of  desorption 
from  the  film  suggest  an  additional  20-30%  enhancement  over 
the  absorption  rale  computed  from  correlations  of  Emmert  and 
Pigford,  implying  that  waves  do  more  than  stretch  and  contract 
the  velocity  profile.  By  imposing  external  disturbances  and 
generating  standing  waves  of  various  amplitudes  and  frequen¬ 
cies  on  stationary  horizontal  films,  Goren  and  Mani  (1968) 
measured  mass  transfer  rates  greater  than  ten  times  the  values 
expected  for  smooth  films,  with  increases  scaling  roughly 
linearly  with  wave  amplitude. 

Data  from  a  variety  of  investigators  published  over  a  three- 
decade  period  were  correlated  to  film  Reynolds  number  and 
Schmidt  number  by  Henstock  and  Hanratty  (1979).  The  data 
reflect  a  variety  of  columns  sizes,  with  Schmidt  numbers  ranging 
from  250  to  1,200  and  Reynolds  numbers  from  100  to  10,000, 
encompassing  laminar  and  turbulent  films.  The  ratio  of  the  mass 
transfer  coefficient  from  their  correlation  to  that  of  the  short 
contact  time  solution,  k is: 


e  -  ~~  «  o.oi  1 1  Ko.707  Jfcy 

*CSCT  V  He 

+  (0.0310 f"-Vife)s}0J0  (6) 

where  kL  is  the  liquid  side  mass  transfer  coefficient  averaged 
over  the  film  length  L  and  hN  is  the  Nussclt  or  time-averaged 
film  thickness.  Equation  4  can  be  used  to  show  that 
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The  enhancement,  E,  depends  only  on  film  Reynolds  number 
and  length,  varying  from  1.15  for  L/h#  of  3,000  and  Reynolds 
number  of  100  to  2.9  for  L/hN  of  10,000  and  Reynolds  number 
of  1 ,000.  An  enhancement  of  nearly  fifteen  times  is  predicted  for 
a  turbulent  film  with  a  Reynolds  number  of  10,000  and  length 
L/hN  of  10,000.  The  additional  interfaciai  area  due  to  the  waves 
over  that  of  a  fiat  interface  has  been  shown  to  be  negligible  over 
a  wide  range  of  flow  rates  (Portalski  and  Clegg,  1971),  so  the 
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enhancement  in  Eq  6  is  not  influenced  by  differences  in 
interfacial  area  due  to  the  presence  of  waves. 

The  modeling  of  gas  absorption  in  falling  films  has  proceeded 
along  four  paths. 

1.  Levich  (1962)  extended  the  short  contact  time  solution  to 
allow  the  surface  velocity  to  vary  with  film  thickness.  Treating 
the  interface  as  a  periodic  array  of  small-amplitude  traveling 
capillary  waves,  Levich  used  the  model  of  Kapitza  and  Kapitza 
(1949)  to  predict  the  surface  velocity  and  concluded  that  the 
enhancement,  £,  would  be  about  1.15.  Ruckenstein  and  Ber- 
bente  (1965.  1968)  also  used  the  Kapitza  hydrodynamic  model 
and  approximated  the  interface  as  small-amplitude  traveling 
waves  described  by  two  Fourier  modes.  The  concentration 
profile  in  the  film  was  approximated  by  a  power  series  in  the 
coordinate  normal  to  the  interface  and  an  enhancement  of  about 
30%  was  predicted.  Since  the  Kapitza  wave  model  describes 
only  capillary  waves  the  results  are  applicable  only  for  very  low 
Reynolds  numbers. 

More  recently,  Barrdahi  (1988)  solved  the  convective- 
diffusion  equation  for  transport  into  wavy  films  at  very  low 
Reynolds  numbers,  near  the  condition  of  wave  initiation.  He 
showed  that  mass  transfer  enhancement  due  to  small  waves 
scales  with  Pe'n.  Howard  and  Lightfoot  (1968)  arrived  at  the 
same  conclusion  by  treating  the  interface  as  a  periodically 
stretching  surface.  Javdani  (1974)  suggested  a  model  for 
wave-induced  concentration  fluctuations  similar  to  simple  eddy 
viscosity  models.  Using  the  Kapitza  model  for  the  velocity 
profile  he  proved  only  that  enhancement  scales  as  Pe'11. 

2.  The  models  discussed  above  are  deficient  in  their  neglect  of 
the  presence  of  large  roll  waves  and  the  substantial  velocities 
normal  to  the  interface  that  can  be  expected  to  accompany 
them.  Attempts  to  relate  the  increased  mass  transfer  to  large 
wave  properties  using  surface  renewal  models  are  summarized 
by  Davies  (1972)  and  Dukler  (1977).  Banerjee  et  al.  (1967) 
proposed  a  renewal  model  which  assumed  that  the  large  waves 
mixed  with  the  substrate  over  which  they  passed,  bringing  fresh 
solution  to  the  interface.  Closing  the  problem  required  a  relation 
between  the  Reynolds  number  and  large-wave  frequency;  this 
was  derived  from  linear  stability  considerations.  Since  linear 
stability  analyses  characterize  only  small  capillary  waves  (Ben¬ 
jamin,  1957),  this  choice  of  relationships  is  questionable.  The 
model  was  modified  by  Brumfield  et  al.  (1975),  incorporating 
new  data  on  large-wave  frequencies  (Telles  and  Dukler,  1970). 

3.  Wave-induced  turbulence  was  suggested  by  Suzuki  et  al. 
(1983)  as  a  means  of  explaining  the  enhanced  mass  transfer 
rates.  The  method  proposed  a  turbulent  diffusivity  for  both  large 
and  small  waves  whose  definition  was  empirically  related  to  the 
size  and  velocity  of  the  waves.  Poor  agreement  between  pre¬ 
dicted  and  measured  transport  rates  was  reported  for  instances 
in  which  this  model  was  applied  to  flows  with  many  large  waves. 

4.  Films  sheared  by  the  surrounding  gas  have  been  studied 
experimentally  .  “icCready  and  Hanratty,  1985)  and  analyti¬ 
cally  (McCready  ct  al.,  1986;  Back  and  McCready,  1988).  It  is 
suggested  that  shear  stress  variations  due  to  gas  flow  around 
waves  induce  normal  velocities  near  the  interface  which  influ¬ 
ence  mass  transfer.  Of  course,  this  mechanism  is  not  applicable 
to  the  case  of  free  falling  films. 

The  experimental  evidence  accumulated  over  the  past  three 
decades  shows  that  mass  transfer  enhancement  due  to  waves  on 
films  cannot  be  explained  with  models  using  any  variant  of  the 
Kapitza  capillary  wave  velocity  profile.  Models  based  on  largc- 


wave-induccd  surface  renewal  suffer  lor  lack  of  a  complete 
characterization  of  wave  structure,  frequency,  and  amplitude 
distribution.  Analyses  based  on  more  robust  hydrodynamic 
models  appear  to  be  the  logical  next  step  toward  reconciling 
experimental  measurements  and  predicted  transfer  rates 

Mass  transfer  from  the  wall 

Few  experimental  studies  of  this  problem  have  been  reported 
Referring  to  Figure  3,  most  experiments  have  been  conducted 
with  a  contaminant  affixed  to  a  tube  wall  or  flat  plate  over  a 
distance  LA.  A  film,  was  allowed  to  flow  over  the  solute,  with 
careful  attention  to  complete  wetting  of  the  surface  Mass 
transfer  rates  were  determined  either  by  weighing  the  plate 
before  and  after  the  experiments  or  by  measuring  the  outlet  bulk 
concentration  of  the  fluid. 

Stirba  and  Hurt  (1955)  attempted  to  relate  the  increased 
transfer  rate  due  to  the  waves  to  a  universal  eddy  diffusivity. 
Experiments  in  2  m  long  tubes  coated  with  organic  acids  over 
lengths  from  1  to  1.5  m  showed  apparent  diffusivities  ranging 
from  three  to  twenty  times  the  molecular  diffusivity,  depending 
on  the  Reynolds  and  Schmidt  numbers.  Reynolds  numbers  in 
the  study  varied  from  300  to  3,000,  and  Schmidt  numbers  from 
600  to  18,000.  No  general  correlation  of  the  results  was 
presented. 

Oliver  and  Alherinos  (1968)  conducted  experiments  in  an 
inclined  channel  at  angles  of  30°  or  less  with  the  horizontal. 
Re  <  200,  and  length  of  about  0.3  m.  They  found  that  transfer 
rates  were  described  adequately  by  a  short  contact  time,  smooth 
liquid  film  theory  analogous  to  that  of  Pigford  for  gas-liquid 
transfer.  However,  the  conditions  of  the  experiment  were  such 
that  large  waves  were  not  present.  Oliver  and  Atherinos  suggest 
that  the  difference  between  their  result  and  that  of  Stirba  and 
Hirt  is  due  to  the  presence  of  large  waves  in  the  latter 
experiments  while  theirs  only  showed  capillary  wave  motion. 

Mass  transfer  from  a  wall  to  a  liquid  is  represented  by  a 
substantial  body  of  literature  because  of  interest  in  the  use  of 
electrochemical  probes  to  measure  wall  shear  stress  (Hanratty 
and  Campbell,  1983).  However,  these  models  all  assume  the 
existence  of  a  very  thin  concentration  boundary  layer  near  the 
wall  through  which  the  velocity  profile  can  be  assumed  to  be 
linear.  This  condition  is  a  very  specialized  one  and  not  of  general 
applicability.  No  analytical  theories  exist  that  can  be  used  to 
explore  the  effect  of  wave  motion  on  this  wall  to  fluid  transfer. 
Speculation  about  the  effect  of  large  waves  on  transfer  rates  is 
possible  through  an  examination  of  the  velocity  profile  near  the 
wall,  where  the  streamwise  velocity  can  be  approximated  using  a 
Taylor  expansion, 

aU.y.t)  -  u(x.O.t)  +  (— )  y  +  .  .  (8) 

\°y  /u.o.ii 

where  u(x,  0,  t)  is  identically  zero,  and  the  derivative  term 
represents  the  wall  shear  stress.  Using  the  continuity  equation, 
the  normal  velocity  is  seen  to  be  approximated  by 

(9) 

This  simple  scaling  shows  that  mass  transport  enhancement  is 
expected  near  the  front  and  rear  or  targe  waves,  where  wall 


1382 


September  1990  Vol.  36,  No.  9 


AICbE  Journal 


hear  stress  values  are  rapidly  changing  (Wasden  and  Duklcr. 
1989a.  b). 

Computational  Procedure 

The  numerical  study  focused  on  simulating  mass  transfer  in  a 
series  of  seven  experimentally  measured  large  wave  shapes 
chosen  from  film  thickness  traces  obtained  at  a  Reynolds 
number  of  880.  The  wave  shapes  chosen  included  both  isolated 
and  interacting  waves  and  were  representative  of  all  large  waves 
existing  on  the  film.  Simulation  of  mass  transfer  through  either 
the  wavy  interface  or  from  the  wall  combined  three  numerical 
algorithms.  First,  the  u  and  v  velocity  fields  were  computed  for 
each  wave  using  a  finite-difference  numerical  algorithm  based 
on  the  TEACH-T  code  (Gosman  and  Idcriah,  1976)  and 
described  elsewhere  (Wasden  and  Dukler,  1989a,  b).  The  hydro- 
dynamic  simulation  was  performed  under  the  assumption  of 
passive  scalar  transport,  that  is,  the  presence  of  the  diffusing 
species  did  not  affect  the  physical  properties  of  the  fluid. 
Boundary  conditions  required  to  specify  the  concentration 
distribution  in  the  flat  film  surrounding  the  large  waves  were 
then  computed  from  a  numerical  solution  of  Eqs.  2  and  3  or  Eqs. 
2  and  S.  The  final  step  in  the  transport  simulation  was  the 
solution  of  the  '-•■'nvective-diffusion  equation  using  the  velocity 
and  spatial  shape  profiles  determined  in  the  hydrodynamic 
simulations. 


described  by  the  solution  of  the  fiat  film  problem,  Eqs.  2  and  3, 
for  a  given  position  of  the  wave  in  the  column,  L  Mass  transfer 
through  the  wavy  interface  is  computed  as  the  wave  passes  a 
location  L  below  the  feed.  A  locally  variable  wave  velocity, 
E„(z),  is  determined  in  the  course  of  the  hydrodynamic  simula¬ 
tion  and  used  to  approximate  the  unsteady  terms  in  Eq  1, 


The  solution  is  then  generated  for  a  steady  problem  in  a 
coordinate  system  moving  with  the  wave.  The  convective- 
diffusion  equation  for  this  planar  flow  is  then  written 


where  the  coordinates  are  shown  in  Figure  4.  At  the  solid 
boundary  ,  y  0,  a  no-flux  condition.  Eq.  3c,  is  imposed  while 
the  free  surface  streamline  corresponds  to  a  saturated  solution, 
Eq.  3b.  The  inlet  condition,  at  r  •»  L,  is  determined  by  a 
numerical  solution  of  Eqs.  2  and  3  for  given  Peclet  number, 
substrate  thickness,  and  fluid  properties.  The  outlet  condition,  at 
r  =  L  -A,  corresponds  to  negligible  sieamwise  concentration 
gradients,  or 


Gas  absorption  through  the  wavy  interface 
The  gas  absorption  problem  simulated  in  this  study  is  illus¬ 
trated  in  Figure  4.  Each  large  wave  is  modeled  as  being 
surrounded  by  a  flat  film  in  which  the  concentration  field  is 


Figure  4.  Absorption  through  large  waves. 


where  X  is  the  wavelength  of  the  large  waves  and  includes  only 
the  sloped  portions  of  the  wave. 

The  solution  of  the  transport  problem  was  generated  from  a 
modified  version  of  the  finite-difference  algorithm  developed  for 
the  hydrodynamic  problem  (Wasden  and  Dukler,  1989b).  In 
order  to  increase  accuracy,  a  novel  version  of  a  quadratic  upwind 
differencing  technique,  QUICK  (Leonard,  1979),  called  NU- 
QUICK-ER  was  developed  for  approximating  convective  terms 
in  both  streamwise  and  normal  directions.  A  mass  source 
calculated  using  the  concentration  gradient  at  the  surface  was 
added  to  the  mass  conservation  equation  for  those  control 
volumes  bordering  the  free  interface,  while  no  mass  was  allowed 
to  leave  the  control  volume  through  the  wall.  Details  of  the 
numerical  method  and  surface  treatment  techniques  appear 
elsewhere  (Wasden,  1989). 

The  concentration  field  within  each  wave  determined  in  this 
way  was  used  to  compute  values  of  the  local  mass  flux  at  the 
interface  through  Fick’s  law.  Of  particular  interest  is  the  mass 
flux  integrated  over  the  wavelength.  The  local  flux  given  by  the 
short  contact  theory  can  be  integrated  over  the  wavelength  to 
yield 


Grid  refinement  studies  showed  that  the  grid  mesh  necessary 
to  fully  resolve  hydrodynamic  details  in  the  large  waves  provided 
sufficient  detail  for  transport  modeling.  Approximately  2,000 
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mesh  cells  were  used  in  each  case.  The  iterative  solution  of  the 
transport  problem  was  continued  until  the  sum  of  mass  residuals 
throughout  the  wave  was  less  than  one  part  in  10*  of  the  total 
mass  transferred  across  the  interlace.  Execution  tune  for  the 
code  was  between  200  and  300  CPU  seconds  on  an  NAS  9000 
mainframe  computer. 

Mass  tansfer  from  the  wall 

Mass  transfer  from  the  wall  was  computed  by  tracking  the 
evolution  of  concentration  profiles  with  time  as  the  waves 
traversed  the  soluble  surface.  The  total  amount  of  mass  accumu¬ 
lated  by  a  wave  was  compared  to  that  numerically  computed  for 
a  smooth  film  having  the  same  time  of  exposure  to  determine  the 
enhancement  due  to  the  presence  of  the  wave. 

This  simulation  procedure  required  an  algorithm  designed 
specifically  for  unsteady  simulations.  The  physical  situation 
described  by  the  model  is  shown  in  Figure  5.  For  this  problem, 
the  unsteady  convective-difTusion  equation  is  given  by  Eq.  1  with 
the  coordinates  as  defined  in  Figure  5.  Time  begins  when  the 
wave  front  reaches  the  upper  end  of  the  active  surface,  where  the 
film  is  assumed  to  be  composed  of  pure  solvent,  Eq.  3a.  Along 
the  soluble  surface,  the  concentration  is  equal  to  the  value  of  a 
saturated  solution,  Eq.  5a,  while  no  flux  is  allowed  at  the  wall  in 
regions  outside  the  active  region,  requiring 

=0,jr<0.x>Z.„  (14) 


Figure  5.  Mass  transfer  from  solid  boundary  into  a  pass¬ 
ing  wave. 


No  mass  is  allowed  to  pass  through  the  wavy  interface,  or 


(15) 


The  problem  formulation  includes  two  parameters,  the  length  of 
the  active  surface,  and  the  Schmidt  number.  As  in  the  case 
of  mass  transfer  to  the  interface,  the  simulation  used  the  same 
rectangular  grid  mesh  as  the  hydrodynamic  simulation;  roughly 
2,000  grid  points  were  specified  for  each  large  wave.  The  small 
grid  spacing  used  near  the  wall  to  resolve  the  wall  shear  stress  in 
the  hydrodynamic  problem  was  suitable  for  the  mass  transfer 
problem  where  the  concentration  varied  rapidly  near  the  wall. 

The  unsteady  term  in  the  convection-diffusion  equation  was 
approximated  with  a  first-order  implicit  finite-difference  tech¬ 
nique  as  recommended  by  Patankar  (1980).  The  unsteady  term 
adds  a  source  term  to  the  discretized  equations  that  is  easily 
assimilated  into  the  control-volume-based  method  on  which  the 
algorithm  was  based.  The  implicit  method  is  unconditionally 
stable  for  the  present  problem,  insuring  that  the  temporal 
behavior  of  the  results  are  reasonable  approximations  to  the  true 
solution.  To  reduce  errors  in  the  solution,  time  step  sizes  of  less 
than  1  ms  were  used. 

The  solution  procedure  began  with  the  concentration  field  set 
to  zero.  As  time  advanced  one  step,  the  wave  moved  a  distance 
determined  by  the  velocity  near  the  wall  for  the  first  grid  cell. 
The  first  streamwise  wall  cell  was  specified  as  having  a  constant 
concentration  boundary  condition,  therefore  allowing  a  flux  of 
solute,  while  all  other  wall  ceils  were  bounded  by  the  solid  wall, 
with  a  no-flux  condition  imposed.  At  this  location  and  time,  the 
concentration  field  was  determined  by  solving  the  system  of 
linear  equations  representing  the  finite-difference  form  of  Eq.  1 
using  an  alternating  direction  implicit  (ADI)  technique.  This 
iterative  technique  was  continued  until  the  solution  had  con¬ 
verged  to  within  one  part  in  10'.  The  wave  was  then  allowed  to 
move  downward  another  time  step.  The  distance  moved  with 
each  time  step  depended  on  the  velocity  associated  with  each 
grid  near  the  wall;  the  algorithm  adjusted  the  time  step  such 
that  each  step  corresponded  to  moving  the  wave  sufficiently  that 
one  more  whole  grid  cell  was  subject  to  a  flux  condition.  The  grid 
mesh  spacing  in  the  streamwise  direction  was  specified  such  that 
smaller  grids  were  assigned  to  regions  in  which  the  velocity  fields 
were  changing  rapidly.  This  grid  spacing  scheme  insured  that 
smaller  time  steps  were  taken  in  these  regions,  and  that  the  mass 
transport  was  accurately  approximated.  Implicit  in  the  simula¬ 
tion  is  the  condition  that  each  wave  continues  to  evolve  in  the 
same  way  that  was  determined  for  a  single  location.  This 
assumption  is  justified  for  an  isolated  wave  changing  very  little 
with  distance.  In  contrast,  evolving  and  interacting  waves  arc 
probably  not  well  described  by  this  criteria  over  a  long  distance. 
The  present  simulation  will  be  limited  to  reasonably  short  active 
surface  lengths  in  order  to  minimize  this  error. 

Concentration  fields  were  determined  for  each  wave  for  fixed 
values  of  the  active  surface  length  and  various  Schmidt  num¬ 
bers.  The  algorithm  required  about  1  CPU  second  (on  an  NAS 
9000  computer)  for  each  time  step,  and  for  most  cases,  200-300 
CPU  seconds  were  sufficient  to  complete  each  simulation. 


Average  Wave  Structure  of  the  Falling  Film 

An  average  wave  structure  was  defined  from  time  records  of 
film  thickness  measurements  taken  at  Re  ~  880  to  determine 


1384 


September  1 990  Vol.  36,  No.  9 


AiChE  Journal 


whether  the  statistically  defined  average  wave,  if  distributed 
along  the  surface,  would  yield  computed  mass  transfer  rates 
comparable  to  experimental  values.  The  film  is  modeled  as 
periodic  average  waves  separated  by  a  rippled  substrate. 

Film  thickness  data  from  which  the  large  wave  profiles  were 
taken  were  examined  to  determine  the  characteristics  of  an 
average  large  wave.  Figure  6  illustrates  the  components  of  the 
proposed  film  structure,  a  series  of  triangular  waves  traveling 
over  a  flat  substrate.  Data  consisting  of  45  s  of  film  thickness 
signals  taken  at  two  locations  and  sampled  at  l  kHz  per  channel 
were  analyzed  using  traditional  statistical  techniques  to  arrive 
at  the  dimensions  and  velocity  of  the  wave  structure. 

The  wave  velocity  was  determined  by  dividing  the  distance 
separating  two  film  thickness  probes  (63  mm)  by  the  time  delay 
of  the  maximum  of  the  cross-covariance  of  the  signals.  The  time 
of  the  first  zero  of  the  film  thickness  autocorrelation  provided  an 
average  wave  frequency,  /,  from  which  the  total  length  of  the 
wave  in  the  time  domain  was  determined  as 


1 

j.  —  tf  +  it,  +  t,  (16) 


where  tf ,  tb,  and  t,  are  the  lengths  (in  a  time  domain)  of  the  front 
of  the  wave,  the  tail,  and  the  separating  substrate.  Averaging 
over  the  ensemble  of  all  measured  waves  separated  by  a 
substrate  thinner  than  the  mean  thickness  provided  average 
values  for  t,  and  tb.  The  extent  of  the  substrate,  is  then 
determined  using  Eq.  16. 

The  average  substrate  thickness,  A,,  was  chosen  as  the  value 
below  which  5%  of  the  film  thickness  points  fell.  This  value 
approximates  the  average  of  the  minimum  values  of  film 
thickness  between  large  waves.  With  this  value,  the  problem  was 
closed  using  the  volume  conservation  equation 


hy 


(hp  -  Ut  +  lf) 

\  2  /  (f,  +  //  +  tb) 


(17) 


which  can  be  solved  for  the  peak  thickness. 


hp  — 


2  (hN  —  h ,)  (/*  +  if  +  t, 


(f»  +  tf) 


+  h. 


(18) 


The  experimental  data  analyzed  with  this  method  had  an 
average  thickness,  hN ,  of  0.365  mm,  an  average  substrate 


thickness.  h„  of  0  260  mm,  and  an  average  peak  thickness.  hr 
equal  to  0.614  mm.  The  peak  to  substrate  ratio  ts  2.36  An 
average  frequency  of  6.9  Hz  was  computed,  along  with  average 
i,of  77.3  ms  and  tb  of  51.0  ms  The  substrate  duration,  f,,  was 
computed  to  be  60,7  ms.  The  wave  velocity  for  the  film  was 
determined  to  be  1.15  m/s.  Streamlines  resulting  from  the 
hydrodynamic  simulation  of  the  average  wave  are  shown  in 
Figure  7. 

Results 

Gas  absorption  through  the  wavy  interface 

Simulations  were  carried  out  for  seven  large  waves  whose 
profiles  were  measured  3.1  m  below  the  feed  and  which  were 
representative  of  all  isolated  and  interacting  waves.  For  each 
wave,  the  Schmidt  number  was  varied  from  250  to  1,000,  while 
the  distance  below  the  feed  ranged  from  1,000  to  6,000  Nusselt 
thicknesses.  This  parameter  space  encompasses  nearly  all  indus¬ 
trially  important  gas-liquid  diffusion  systems  and  column 
heights.  For  each  of  the  parameter  pairs,  a  concentration  field 
was  determined.  The  local  variation  of  mass  flux  could  then  be 
computed  from  the  gradient  of  this  concentration  at  the  inter¬ 
face  and  the  total  flux  through  the  interface  determined  from 
the  integration  of  this  local  flux  along  the  wave.  Two  waves  will 
be  discussed  in  detail:  an  isolated  wave  with  peak/substrate 
thickness  of  approximately  three  and  an  interacting  wave. 
Streamline  maps  for  these  waves  are  given  in  Figures  8  and  9. 

Concentration  fields  computed  for  the  isolated  large  wave  are 
shown  in  Figure  10  for  Schmidt  numbers  of  500  and  1,000  at  a 
distance  of  4,500  mean  film  thicknesses  below  the  feed  (1 .63  m 
for  a  Reynolds  number  of  880).  These  conditions  are  compara¬ 
ble  to  those  investigated  experimentally  by  Emmert  and  Pigford 
( 1954) .  The  concentration  profiles  are  presented  as  contour  plots 
with  intervals  of  0.1  saturated  concentration  units  dividing 
successive  contour  lines.  Both  concentration  profiles  show  signif¬ 
icant  deviations  from  parallel  isoconcentration  lines  predicted 
by  the  short  contact  time  theory  (SCTT).  Normal  (v)  velocities 
within  the  wave  peak  carry  solute  from  the  surface  into  the 
wave.  Near  the  front  and  rear  of  the  wave,  normal  velocities 
resulting  from  the  wave  peak  interacting  with  the  substrate 
force  the  solute  into  the  substrate.  The  effect  of  the  Schmidt 
number  is  seen  in  the  extent  to  which  the  concentration  field 
penetrates  the  wave  substrate;  for  a  higher  Schmidt  number, 
this  penetration  occurs  to  a  lesser  extent. 


or 
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Figure  6.  Average  film  structure  dimensions. 


Figure  7.  Streamline  map  for  average  large  wave. 
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Figure  8.  Streamline  map  for  representative  Isolated 
wave. 


A  profile  of  the  local  mass  flux  {for  a  saturated  concentration 
of  1  kg/mJ)  for  the  isolated  wave  is  presented  in  Figure  11  for  a 
Schmidt  number  of  500.  The  absence  of  normal  velocities  near 
the  interface  greatly  affects  the  local  flux  near  the  front  and  rear 
stagnation  points,  where  the  local  fluxes  approach  those  ex¬ 
pected  for  a  purely  diffusive  situation.  In  contrast,  near  the  front 
and  behind  the  peak  of  the  wave,  convection  increase  the  local 
flux  by  several  times  over  the  diffusive  limit.  Comparison  of 
Figures  10  and  II  shows  that  the  maxima  in  local  flux  values 
correspond  to  regions  in  which  the  wave  peak  interacts  with  the 
slowly  moving  substrate.  This  hydrodynamic  process  is  also 
responsible  for  the  dramatic  increase  in  wall  shear  stress  over 
that  for  parallel  flow.  Concentration  and  local  flux  profiles  for 
the  average  wave  were  similar  to  those  computed  for  the  isolated 
wave.  Previous  models  that  suggest  surface  renewal  as  a  result  of 
large  waves  interacting  with  the  substrate  are  now  seen  as 
clearly  inconsistent  with  the  computed  concentration  profiles. 

Figure  12  shows  the  concentration  profile  in  an  interacting 
wave  for  a  Schmidt  number  of  500  appearing  4,500  mean  film 
thicknesses  below  the  feed.  The  local  flux  profile  for  this 
interacting  wave  is  shown  in  Figure  13.  The  profile  appears  to 
approximate  the  sum  of  the  profiles  obtained  from  isolated 
waves.  The  extra  stagnation  points  create  additional  extrema  in 
the  profile,  with  the  minimum  occurring  in  the  trough  between 
the  waves  as  a  result  of  the  stagnation  region  separating  the  two 
counterrotating  recirculation  zones.  The  presence  of  the  Stag- 


Distance.  mm 

Figure  10.  Concentration  profiles  in  isolated  wave  4,500 
mean  Mm  thicknesses  below  feed  for  Schmidt 
nos.  500  and  1,000. 


nant  zone  separating  the  peaks  dampens  the  effect  of  the 
interacting  waves  on  mass  transport.  Compared  to  the  SCTT 
prediction,  enhancements  arc  close  to  those  associated  with 
isolated  waves. 

The  integral  of  the  local  flux  along  the  wave  was  normalized 
by  that  expected  for  SCTT  as  given  in  Eq.  13.  The  predicted 
mass  transfer  to  an  isolated  wave  was  then  computed  by 
weighting  the  flux  through  the  wavy  interface  and  the  rippled 
substrate.  Since  the  substrate  was  covered  with  capillary  waves 
it  was  modeled  using  an  enhancement  of  30%  as  given  by  the 
model  of  Ruckenstein  and  Berbente  (1968).  The  enhancement 
for  the  wave  and  its  associated  substrate  was  then  computed 
from: 


Total  flux 

'/+'»  I 

Flux  due  to  wr  /e 

SCTT  prediction 

[tf+h  +  t.l 

SCTT  prediction 

30)  (19) 


Enhancement  factors  computed  for  the  seven  waves  whose 
profiles  were  measured  are  presented  in  Figure  14  for  several 
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Figure  12.  Concentration  profiles  in  Interacting  wave 
4,500  mean  film  thicknesses  below  feed  for 
Sc  —  500. 


locations  below  the  entry.  It  should  be  recognized  that  the  wave 
profile  and  velocity  vary  along  the  length  of  the  film  (Zabaras 
and  Dukler,  1988)  while  for  these  computations  the  profile 
measured  at  LjhN  =  8,550  (3.1  m)  was  used  for  the  mass 
transfer  computation  at  all  locations.  The  numerical  analysis 
predicts  mass  transfer  rates  1. 7-3.5  times  greater  than  the  short 
contact  time  predictions.  A  comparison  with  experimental  data 
as  correlated  by  Henstock  and  Hanratty  (1979)  is  included  in 
this  figure.  The  experiments  represent  mass  transfer  rates 
averaged  from  L/hN  of  zero  to  the  indicated  value  while  the 
computation  describes  the  local  rate  as  the  wave  passes  a 
particular  L/h„.  Similar  qualitative  trends  in  enhancement  are 
clearly  seen  in  both  computed  and  measured  values. 

Figure  15  compares  the  enhancement  computed  for  the 
average  wave  uniformly  distributed  across  the  length  of  the 
interface  with  the  experimental  enhancement  given  by  the 
Henstock  and  Hanratty  (1979)  correlation.  The  agreement  in 
magnitude  and  trend  is  now  semiquantitative.  Furthermore,  the 
discrepancies  can  readily  be  understood.  For  short  columns  the 
waves  are  smaller  and  the  rate  of  mass  transfer  would  be  less 
than  that  computed  for  the  average  film  structure  in  the  figure. 
At  large  values  of  LjhN  the  existence  of  waves  above  this 
location  changes  the  inlet  condition  used  in  the  computations  to 
one  of  a  more  uniform  concentration  as  a  result  of  the  successive 


Figure  13.  Local  mass  flux  values  for  interacting  wave  for 
conditions  of  Figure  12. 


Figure  14.  Computed  and  observed  mass  transfer  en¬ 
hancement  for  transfer  through  wavy  inter¬ 
face  at  Sc  =  500. 


mixing  induced  by  the  upstream  waves.  This  would  lead  to 
larger  concentration  gradients  at  the  surface  and  consequently  a 
higher  mass  transfer  rate  than  predicted  by  the  simulation. 
Inclusion  of  these  effects  would  produce  even  closer  agreement 
between  the  model  and  experiment.  Once  the  ability  to  predict 
the  average  profile  as  a  function  of  position  along  the  film  is 
developed  these  factors  can  readily  be  incorporated  into  the 
computation. 

Mass  transfer  from  the  wall 
Simulations  of  the  mass  transfer  in  the  two  waves  discussed  in 
the  previous  section  will  be  presented.  Concentration  fields, 
average  concentration  values  [across  the  film,  from  y  —  0  to 
y  =  h(x, »)],  and  total  mass  accumulation  were  determined  for 
fixed  values  of  the  active  surface  length  and  various  Schmidt 
numbers.  Mass  transfer  enhancement  due  to  large  waves  was 
determined  by  computing  the  total  mass  within  the  large  wave 
after  it  had  passed  over  the  active  surface  and  comparing  it  to 
the  amount  of  mass  accumulated  by  a  flat  film  of  thickness 
which  had  been  in  contact  with  the  active  surface  for  the  same 
amount  of  time.  The  active  surface  length  for  the  present  results 
was  800  average  film  thicknesses,  or  roughly  0.29  m  for  a 
Reynolds  number  of  880.  This  value  represents  about  one 
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Figure  15.  Mas*  transfer  enhancement  computed  for  av¬ 
erage  film  structure  compared  to  observed 
values  for  Sc  =  500. 
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wavelength.  Results  for  Schmidt  numbers  of  200  and  1,200  are 
presented.  Fluid  physical  properties  were  taken  as  identical  to 
those  used  for  the  hydrodynamic  simulation,  leading  to  Peclct 
numbers  ranging  from  1.76  x  105  to  1.056  x  10*. 

Average  concentration  profiles  (normalized  by  the  saturated 
wall  concentration)  for  the  isolated  large  wave  are  presented  in 
Figure  16.  The  increase  in  average  concentration  ahead  of  the 
peak  is  due  to  the  strong  normal  velocities  at  the  wall  in  this 
region,  which  result  from  the  rapid  change  in  wall  shear  stress. 
The  average  concentration  profiles  appear  to  be  affected  most  by 
the  strong  hydrodynamic  interactions  between  the  wave  peak 
and  body  with  the  surrounding  fiat  film.  Other  hydrodynamic 
processes  caused  by  the  wavy  interface  apparently  do  not 
penetrate  to  a  level  sufficient  to  alter  the  transport  processes  at 
the  wall. 

Further  evidence  of  the  strong  normal  velocities  near  the  wall 
is  found  in  the  profiles  of  near  wall  concentrations,  Figure  16b. 
The  near  wall  concentrations  are  those  that  are  found  in  the  first 
cell  inside  the  boundary  and  are  proportional  to  the  concentra¬ 
tion  gradient  at  the  wail.  These  values  are  significantly  lower 
directly  under  the  steepest  part  of  the  wave  front,  which 
corresponds  to  the  region  in  which  the  shear  stress  is  a 
maximum.  Average  and  near  wall  concentrations  are  closest  in 
magnitude  near  the  front  and  under  the  peak  of  the  wave.  In 
contrast,  the  surrounding  regions  of  substrate  are  characterized 
by  near  wall  and  average  concentrations  that  differ  by  as  much 
as  a  factor  of  ten. 

Increasing  the  Schmidt  number  reduces  mass  diffusivity. 


increasing  the  resistance  to  mass  transport  into  the  waves. 
However,  the  presence  of  significant  normal  velocities  near  the 
wall  reduces  this  effect.  For  large  Peclct  numbers,  increasing  the 
Schmidt  number  sixfold  (from  200  to  1 .200)  would  be  expected 
to  decrease  the  average  concentration  by  the  same  factor  if  no 
normal  velocities  were  present.  However,  Figure  16a  shows  that 
the  average  concentration  is  reduced  by  no  more  than  three 
when  the  hydrodynamic  effects  are  included.  In  a  flat  film,  the 
near  wall  concentrations  would  be  expected  to  scale  inversely 
with  Schmidt  number,  as  the  flux  is  directly  proportional  to  the 
diffusivity.  Figure  16b  shows  that  the  near  wall  concentration 
decreases  by  no  more  than  one-half  when  the  Schmidt  number  is 
increased  sixfold,  further  illustrating  the  extent  to  which  normal 
velocities  alter  the  mass  transfer  resistance. 

Flow  within  the  interacting  wave  produces  the  concentration 
profiles  shown  in  Figure  17.  The  presence  of  two  interacting  but 
separate  closed  recirculation  regions  within  the  peaks  signifi¬ 
cantly  increases  mass  transfer  through  the  wavy  surface  and 
apparently  has  the  same  effect  on  transport  from  the  wall. 
Pronounced  changes  in  the  average  concentration  accompany 
increasing  mass  diffusivity  in  the  interacting  wave.  While  the 
near  wall  concentration  variation  is  magnified  by  low  diffusivity, 
the  variation  in  average  values  suggests  that  once  the  mass  is 
transported  into  the  peak  regions,  diffusion  becomes  important, 
as  the  average  concentration  profile  associated  with  a  lower 
diffusivity  is  smoother  than  expected. 

Mass  transfer  enhancements  for  three  interacting  waves,  the 
isolated  wave,  and  the  average  film  are  presented  in  Figure  18. 
Experimentally  determined  enhancements  from  Stirba  and 


Figure  16.  Variation  of  concentration  with  Aim  thickness  Figure  17.  Variation  of  concentration  with  Aim  thickness 
and  Schmidt  number  for  isolated  wave.  and  Schmidt  number  for  Interacting  wave. 

(»)  Average  concentration  (a)  Average  concentration 

(b)  Rear  wall  concentration.  (b)  Near  wall  concentration. 
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Figure  18.  Computed  and  observed  mass  transfer  en¬ 
hancement  for  transfer  from  wall  for  Re  = 
880. 


Hurt  (1955)  at  Reynolds  numbers  between  850  and  900  are 
included  for  comparison.  Values  for  the  average  film  fall 
between  the  large  waves  and  agree  reasonably  well  with  the 
data  at  higher  Schmidt  numbers.  At  lower  Schmidt  numbers, 
the  simulations  appear  to  overprcdict  the  enhancement.  The 
discrepancy  results  from  two  experimental  limitations.  Stirba 
and  Hurt  report  difficulties  in  fully  wetting  the  surface  of  the 
benzoic  acid  with  water  (Schmidt  number  609),  which  would 
decrease  the  accumulation  by  the  wavy  film.  For  alcohol- 
organic  acid  systems  (Schmidt  numbers  >  800)  no  wettability 
problems  were  encountered  and  the  comparison  is  correspond¬ 
ingly  more  favorable. 

In  the  limit  of  very  small  Schmidt  numbers,  the  enhancement 
due  to  the  complex  hydrodynamics  vanishes,  but  for  Schmidt 
numbers  as  low  as  200  this  enhancement  may  still  be  as  large  as 
500%.  As  the  Schmidt  number  increases  above  roughly  1,200 
the  enhancement  associated  with  large  waves  begins  to  decrease 
due  to  the  inability  of  the  fluid  to  diffuse  into  the  regions  of  the 
flow  in  which  normal  ( v )  velocities  exist.  It  is  expected  that  for 
Schmidt  numbers  of  industrial  interest,  normally  less  than 
5,000,  the  enhancement  due  to  large  wave  hydrodynamics  will 
be  at  least  500%. 

Conclusions 

The  interface  of  a  falling  liquid  film  is  covered  with  a  random 
array  of  small  and  large  waves  interacting  in  a  complicated 
manner,  including  many  with  peak  thicknesses  several  times  the 
mean.  Neglecting  the  influence  of  the  large  waves  causes  most 
models  to  seriously  underpredict  experimentally  observed  trans¬ 
fer  rates.  Numerical  simulations  of  the  hydrodynamics  within 
large  waves  presented  previously  (Wasden  and  Dukler  1989a,  b) 
exposed  mechanisms  through  which  transport  would  be  aug¬ 
mented.  Hydrodynamic  data  obtained  from  this  procedure  were 
used  in  the  simulation  of  passive  transport  in  these  large  waves. 
The  importance  of  the  complex  hydrodynamics  within  the  waves 
has  been  demonstrated  by  the  semiquantitative  agreement 
between  computed  and  measured  transfer  rates  at  both  inter¬ 
faces.  Future  modeling  of  transport  processes  in  wavy  films  must 
include  provisions  for  predicting  the  significant  normal  velocities 
at  each  interface  and  relating  them  to  the  interfacial  structure. 

The  random  interface  has  been  approximated  by  an  interface 
composed  of  periodically  occurring,  statistically  average  waves 


interspersed  among  smaller  ripples.  Mass  transfer  rates  com¬ 
puted  for  such  an  array  are  1.5  to  2.5  times  greater  than  that 
predicted  for  gas  absorption  into  nonwavy  films  and  four  to 
seven  times  greater  than  the  flat  film  prediction  for  a  dissolving 
wall.  These  enhancements  are  in  reasonable  agreement  with 
experimental  observations. 
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Notation 

c  —  local  contaminant  concentration,  kg/m1 
D  —  gas-liquid  or  liquid-liquid  diffusion  coefficient,  m’/s 
E  —  mass  transfer  rate  enhancement,  Eq.  6 
/  —  average  wave  frequency,  Hz 
g  —  acceleration  of  gravity,  m/s1 
A  =  local  film  thickness,  mm 
A,  —  substrate  thickness  of  average  wave,  mm 
A,  =  peak  thickness  of  average  wave,  mm 
hN  -»  time  average,  or  Nusselt  film  thickness,  mm 
kL  **  liquid-side  mass  transfer  coefficient,  m/s 
L  =  length  over  which  gas  contacts  a  falling  film,  m 
La  =  length  of  contaminated  portion  of  wall,  m 
n  =  coordinate  normal  to  wavy  interface,  m 
Nscn  ~  mass  transfer  rate  predicted  by  short  contact  time  theory, 
kg/m  •  s 

Pe  =  Peclet  number.  Re  ■  Sc 

Q  —  liquid  film  flow  rate  per  unit  perimeter,  f/A  A„,  mJ/s 
Re  —  film  Reynolds  number,  4g/» 

Sc  —  Schmidt  number,  e/D 
lt  =  duration  of  average  wave  back,  ms 
tf  =  duration  of  average  wave  front,  ms 
t,  «*  duration  of  substrate  separating  average  waves,  ms 
u  «.  local  streamwise  velocity,  m/s 
UN  —  time  average,  or  Nusselt  film  velocity,  m/s 
v  ■=  local  velocity  normal  to  boundary,  m/s 
Vm  ■=  wave  velocity,  m/s 
x  =  axial  coordinate  in  lab  frame,  m 
y  —  coordinate  normal  to  boundary,  m 
i  =  axial  coordinate  fixed  on  wave,  m 

Greek  letters 

p  =>  liquid  density,  kg/m1 
v  =  liquid  kinematic  viscosity,  mJ/s 
p  =  liquid  absolute  viscosity,  kg/m  s 
t„  *=  wall  shear  stress,  N/m1 
X  =  length  of  sloped  sections  of  large  wave,  m 
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Insights  into  the  Hydrodynamics 
of  Free  Failing  Wavy  Films 


Three  isolated  waves  of  differing  amplitude  and  shape  were  selected 
from  experimental  measurements  of  a  falling  liquid  film  at  Re  -  880  for 
study  using  an  algorithm  developed  for  solution  of  the  Navier-Stokes 
equations.  The  method  computes  the  velocity  and  pressure  fields  as 
well  as  the  velocity  of  the  wave.  The  results  show  that  large  streamwise 
accelerations  exist  along  with  regions  of  recirculating  flow  in  a  moving 
coordinate  system.  These  features  can  explain  the  enhanced  rates  of 
heat  and  mass  transfer  observed  in  wavy  film  flow.  Computed  wave 
velocities  and  wall  shear  stress  were  in  reasonably  good  agreement 
with  measurements.  Wave  velocity  is  shown  to  be  sensitive  to  small 
variations  in  the  wave  shape  and  explains  the  apparent  random  varia¬ 
tion  of  wave  velocity  with  amplitude  that  has  been  observed  experimen¬ 
tally.  This  numerical  experiment  points  to  the  shortcomings  of  the  many 
methods  used  to  model  large  waves  on  falling  films  that  have  been 
based  on  parabolic  velocity  profiles. 
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Introduction 

Thin  liquid  films  falling  under  the  influence  of  gravity  along 
solid  surfaces  are  encountered  in  a  wide  variety  of  industrial 
process  equipment,  including  wetted-wall  absorbers,  falling- 
film  chemical  reactors,  condensers,  and  vertical  tube  evapora¬ 
tors.  Reliable  design  of  these  processes  depends  on  the  ability  to 
accurately  predict  the  transport  rates  of  heat  and  mass  to  the 
flowing  film.  At  flow-rates  of  industrial  interest,  falling  films 
(even  in  the  absence  of  gas  flow)  evolve  to  a  highly  irregular 
wavy  interface.  Figure  1  displays  a  short  time  trace  of  such  a 
falling  film.  The  surface  is  covered  by  a  complc  \  array  of  large 
and  small  waves  moving  over  a  substrate  that  is  less  than  the 
mean  film  thickness.  The  large  waves,  which  range  in  amplitude 
from  two  to  five  times  the  substrate  thickness,  carry  a  large  frac¬ 
tion  of  the  total  mass  flowing,  and  are  speculated  to  control  the 
rate  of  transport  (Dukler,  1977).  Before  the  heat  or  mass  trans¬ 
fer  rates  to  such  films  can  be  modeled  it  will  be  necessary  to 
understand  the  velocity  distributions  that  exist  within  these 
waves,  as  well  as  their  evolution.  The  present  work  focuses  on 
the  former  question. 

Making  reliable  experimental  measurements  of  the  velocity 
distribution  in  the  films  is  exceedingly  difficult  due  to  the 
extremely  small  film  heights  1  mm),  very  short  passage  time 
of  each  wave  («60  ms)  and  the  random  location  of  the  wave 
height,  as  seen  in  Figure  1 .  Nonintrusive  methods  such  as  LDA 
(Laser  Doppler  Annenometer)  do  not  provide  sufficiently  fine 
resolution  to  investigate  velocity  profiles.  Thus,  experimental 
measurements  appear  limited  to  the  time  variation  of  wall  shear 


stress  and  film  thickness.  As  a  result,  analytical  models  have 
been  developed  in  the  absence  of  hard  data  on  the  true  flow  con¬ 
ditions  that  appear  to  exist  in  the  waves. 

Most  analytical  models  extend  the  concepts  advanced  by 
Kapitza  (1964)  based  on  the  use  of  a  parabolic  velocity  profile 
and  assuming  that  the  streamwise  hydrodynamic  variables  scale 
with  the  wavelength.  In  examining  various  models  developed  to 
that  date,  Dukler  (1972)  concluded  that  all  failed  to  accurately 
represent  any  measured  characteristics  of  the  wave  except  at 
Reynolds  numbers  well  below  those  of  industrial  interest. 

Maron  et  al.  (1985)  treated  isolated  waves  as  a  series  of  seg¬ 
ments,  each  having  a  different  type  of  velocity  distribution 
depending  on  the  physics  of  the  region.  In  the  substrate,  a  para¬ 
bolic  velocity  profile  was  adequate,  while  the  flow  under  the 
front  of  the  wave  was  assumed  to  be  fully  mixed.  The  slowly 
varying  wave  back  was  described  with  a  boundary  layer  model. 
Upon  matching  these  solutions  at  the  segment  boundaries,  it 
was  possible  to  predict  wave  mean  characteristics  (height, 
length,  velocity,  substrate  thickness)  in  reasonable  agreement 
with  the  values  measured  by  Zabaras  (1985).  The  model  was 
fitted  with  a  limited  amount  of  data  from  experimental  mea¬ 
surements  and  it  failed  to  explain  the  large  variation  observed  in 
individual  wave  amplitudes  and  lengths. 

Modeling  the  wavy  film  flow  by  a  direct  solution  of  the  Nav¬ 
ier-Stokes  equations  is  hampered  by  numerical  stiffness  im¬ 
posed  by  the  stress-free  interface;  as  a  result,  convergence  is  dif¬ 
ficult  except  at  the  lowest  flow  rates.  Bach  and  Villadsen  (1984) 
explored  the  application  of  a  finite-element  scheme  to  the 
unsteady  problem  of  waves  developing  from  initial  perturbations 
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Figure  1.  Film  thickness  time  trace. 


on  the  smooth  film  for  Reynolds  numbers  up  to  100.  Their  work 
predicted  that  the  equilibrium  condition  would  consist  of  waves 
having  one  general  shape,  a  condition  contrary  to  experimental 
fact  even  at  film  Reynolds  numbers  as  low  as  1.  The  film  Rey¬ 
nolds  number  is  defined  as  Re  -  4  Q/v,  where  Q  is  the  mass  flow 
rate  per  unit  perimeter,  and  v  is  the  kinematic  fluid  viscosity. 
Kheshgi  and  Scriven  (1987)  applied  a  finite-element  technique 
to  a  problem  with  periodic  boundary  conditions  in  the  flow 
direction,  and  verified  the  evolution  of  infinitesimal  distur¬ 
bances  as  predicted  by  Orr-Sommerfeld  analyses.  Their  work 
was  limited  to  low  flow  rates,  and  failed  to  generate  waveforms 
comparable  to  those  observed  experimentally  for  fully  developed 
flow. 

In  the  absence  of  analytical  models  for  velocity  profiles  that 
appear  to  represent  reality,  and  due  to  the  absence  of  suitable 
experimental  methods  for  measuring  these  profiles,  a  series  of 
numerical  experiments  were  undertaken.  Wave  shapes,  wall 
shear  stress  profiles,  and  wave  velocities  were  measured  in  our 
laboratory  for  a  film  Reynolds  number  of  880,  chosen  to  insure 
significant  inertial  forces  while  remaining  viscous  in  nature.  A 
novel  method  of  solving  free  surface  flows  was  developed,  using 
experimentally  determined  film  thickness  data  for  large,  iso¬ 
lated  waves  to  solve  for  the  position  of  the  free  interface.  The 
results  of  these  computations  demonstrate  the  complex  depen¬ 
dence  of  velocity  distributions  on  wave  shape,  and  represent  an 
early  step  toward  realistic  modeling  of  large  waves. 

Experimental  Procedure 
Flow  loop 

For  fully  developed  wavy  film  flow,  film  thickness  and  wall 
shear  stress  data  were  collected  in  a  50.8  mm  ID  vertical  test 
section  in  a  flow  loop  described  by  Zabaras  el  al.  (1986).  After 
being  pumped  through  a  calibrated  rotameter,  the  aqueous  solu¬ 
tion  entered  the  column  through  an  annulus  whose  inner  wall 
consists  of  a  stainless  steel  porous  sinter  having  100  *im  pore 
size.  Combined  with  careful  leveling  of  the  column  prior  to  data 
collection,  this  entry  section  insured  minimal  deviations  from 
axisymmetric  flow  and  produced  a  smooth  inlet  flow.  The  mea¬ 
suring  station  was  located  3.1  m  below  this  entry  section. 

Measuring  station  and  measurement  techniques 

The  measuring  station,  shown  in  Figure  2,  is  patterned  after 
that  described  by  Zabaras  et  at.  The  removable  section  allows 


simultaneous  measurement  of  film  thickness  and  wall  shear 
stress  at  one  location  and  of  thickness  at  another.  The  station 
was  constructed  of  the  same  material  as  the  flow  loop,  and  was 
carefully  machined  to  insure  a  smooth  transition  to  the  station. 

Film  thickness  probes  consisted  of  twin  parallel  platinum- 
13%  rhodium  wires  of  0.05  mm  dia.,  spaced  2.5  mm  apart, 
which  penetrated  the  flow.  As  described  in  detail  by  Brown  et 
al.,  (1978),  a  linear  relation  exists  between  the  resistance  of  the 
film  between  the  wires  and  the  film  thickness.  Calibration  pro¬ 
ceeded  by  setting  the  measuring  station  horizontal,  blocking  the 
ends,  and  introducing  different  fluid  levels,  determined  to  within 
10  #im  by  using  a  cathetometer,  followed  by  measurement  of  the 
resulting  resistance.  Downstream  electronics  for  converting  this 
resistance  to  a  DC  voltage  signal  are  described  elsewhere  (Za¬ 
baras  et  al.,  1986).  Conductance  of  the  fluid  was  monitored 
closely  at  all  times  during  the  calibration  and  data  collection 
procedures  to  insure  proper  correction  of  any  thermally  induced 
conductance  drift. 

Wall  shear  stress  measurements  were  based  on  the  electro¬ 
chemical  mass  transfer  method  described  by  Hanratty  and 
Campbell  ( 1 983).  For  the  present  series  of  measurements,  the 
iodine/tri-iodide  system  was  chosen.  The  working  solution  con¬ 
tained  0.1M  KI  and  0.004M  l2(s)  in  demineralized  water,  and 
was  replaced  every  2  h  to  minimize  errors  due  to  iodine  evapora¬ 
tion.  A  dry  nitrogen  atmosphere  was  used  in  the  flow  loop  to 
minimize  oxygen  saturation  of  the  solution.  Fluid  properties  ai 
25°C  are:  density,  1,010  kg/m5;  absolute  viscosity,  8.50  x  10  4 
kg/m  ■  s;  and  surface  tension,  7  12  x  10  1  N/m.  The  cathode 
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for  this  system  consists  of  a  flush-mounted  strip  of  platinum  foil, 
0.075  mm  (in  the  flow  direction)  by  1  mm  wide,  embedded  in 
Plexiglas  to  insure  electrical  isolation.  By  measuring  the  current 
produced  by  an  electrochemical  reaction  at  the  surface  of  the 
cathode,  the  wall  shear  stress  at  that  location  is  determined.  For 
the  redox  reaction 

1,"  +  2e— >31  (cathode) 

31  — *  lj  -t-  2e  (anode) 

a  concentration  boundary  layer  develops  on  the  cathode  surface, 
which  is  polarized  at  -0.8  VDC  to  insure  the  concentration 
approaches  zero.  For  the  iodine  system,  the  range  of  polarization 
voltage  is  quite  broad,  insuring  that  large  increases  in  flow  rate 
will  not  deplete  the  electron  source  at  the  cathode.  Details  con¬ 
cerning  the  downstream  electronics  and  calibration  associated 
with  this  measurement  technique  are  found  elsewhere  (Zabaras 
et  al.,  1986). 

It  is  now  recognized  (Mao  and  Hanratty,  1985)  that  the 
response  of  the  electrochemical  probe  is  highly  dependent  on  the 
nature  of  the  “input”  wall  shear  stress.  For  the  ionic  system 
employed  in  this  study,  errors  in  both  phase  and  magnitude  are 
expected  to  be  small  due  to  the  large  ( 101  s~ ')  mean  velocity  gra¬ 
dient,  small  cathode  surface  area,  and  large  Schmidt  number 
(v/D)  of  the  fluid  (=780).  The  relationship  given  by  Hanratty 
and  Campbell  (1982)  between  cathode  current  and  wall  shear 
stress  was  used  in  this  study,  as  the  frequencies  in  the  data  were 
sufficiently  low  to  allow  the  use  of  a  quasi-steady  analysis. 

Data  collection,  processing,  and  analysis 

Voltage  signals  from  two  film  thickness  probes  and  the  wall 
shear  stress  probe  were  first  low-pass  filtered  at  i  kHz,  then  fed 
to  a  microcomputer-based  analog  to  digital  (a/D)  converter. 
Each  signal  was  digitized  at  1  kHz  by  a  Data  Translation  12  bit 
A/D  converter  installed  in  a  DEC  Micro  1 1  /73  microcomputer. 
The  data  set  comprised  1  min  of  data,  and  was  stored  on  the  sys¬ 
tem  Winchester  disk  prior  to  applying  calibration  curves  and 
writing  the  data  to  magnetic  tape  for  further  analysis.  Digitiza¬ 
tion  and  collection  errors  are  expected  to  be  negligible  for  all 
data,  while  calibration  errors  for  the  film  thickness  measure¬ 
ment  are  expected  to  be  less  than  3%.  Errors  inherent  in  apply¬ 
ing  steady  state  wall  shear  stress  calibration  curves  depend  on 
the  nature  of  the  input  signal,  requiring  separate  examination  of 
individual  results.  Zabaras  (1985)  reports  estimated  errors  of 
less  than  7%  for  this  technique. 

Film  thickness  and  wall  shea'-  stress  data  was  examined  to 
locate  isolated  waves,  defined  by  a  wave  having  a  peak  to  sub¬ 
strate  thickness  ratio  greater  than  2,  and  surrounded  by  at  least 
one  wavelength  of  reasonably  flat  film.  For  the  sequences  of  raw 
data,  three  representative  waves  of  various  dimensions  were 
chosen  for  computational  domains.  For  each  case,  a  nominal 
wave  velocity  was  determined  from  the  time  necessary  for  the 
wave  to  travel  from  the  upper  to  lower  film  thickness  probes. 

Numerical  Method 

Solution  cf  frcc-boundary  problems  requires  methods  for 
both  the  solution  of  the  governing  momentum  equations  and 
■ !  »pe  determination.  The  velocity  and  pressure  fields  within  the 
v.  ive  wer<  ..  cmined  by  solving  the  Navier-Stokes  equations  in 
r*  imitivc  >  i*lc  form.  F,jr  a  film  Reynolds  number  of  880,  the 


wave  thickness  generally  was  less  than  1%  of  the  pipe  radius, 
and  therefore  a  two-dimensional  Cartesian  coordinate  system 
was  chosen.  The  transformation  of  time  traces  of  film  thickness 
to  this  coordinate  system  comprised  the  shape  determination 
portion  of  the  overall  algorithm.  The  common  method  of  com¬ 
puting  the  position  of  a  free  interface,  h(x),  is  to  determine  the 
value  of  the  film  thickness,  h,  at  given  values  of  the  streamwise 
variable,  x.  The  present  method  inverts  the  process:  for  given, 
measured  values  of  h,  wc  find  the  values  of  x  that  result  in  h{x) 
satisfying  all  cr  the  free  interface  boundary  conditions.  To 
insure  accurate  representation  of  the  interfacial  pressure,  a 
fourth-order-accurate,  divided-differences  scheme  was  used  to 
compute  the  curvature  of  the  interface: 

h„/(\  +  hl)>'2  (1) 

Initially,  the  waves  were  modeled  as  though  their  shape 
remained  constant  with  time;  these  waves  are  termed  “solitary.” 
The  new  streamwise  coordinate,  z,  is  fixed  on  the  wave,  and 
originated  at  the  front  of  the  wave.  The  film  thickness  profile  in 
the  time  domain,  h(t,),  was  converted  to  the  length  domain, 
fi(z,),  through  the  transformation 

A  -  r.  +  VWD,  -  O  (2) 

for  i  ranging  from  1  to  the  number  of  film  thickness  points  in  the 
isolated  wave.  In  this  manner,  the  wave  profile  was  “stretched  ’ 
for  use  as  a  computational  domain,  and  time  was  removed  from 
the  problem.  For  this  coordinate  system,  the  wave  remains  fixed, 
and  the  wall  moves  upw  ird  at  a  constant  speed  given  by  K„,  the 
wave  velocity  for  the  solitary  wave. 

It  is  useful  to  define  a  new  streamwise  velocity  component, 

«(z,  >•)  -  u\x,  y)  +  Vw,  (3) 

where  u'  (x,  y)  is  the  streamwise  velocity  in  a  coordinate  system 
fixed  on  the  wall.  The  governing  equations  for  this  viscous, 
incompressible,  and  isothermal  flow  relative  to  the  moving  wave 
become 

uu,  +  vuy  -  -PJp  +  eA u  +  g  (4) 

uv,  +  vv,  =  -Py/p  +  kAd  (5) 

u.  +  v f  -  0  (6) 

where  v  is  the  velocity  in  the  normal  (y)  direction,  P  is  the  pres¬ 
sure,  g  represents  gravitational  acceleration,  and  v  and  p  are  the 
kinematic  viscosity  and  density  of  the  fli  id,  respectively.  At  the 
stres;-free  interface,  y  -  A(z),  tangential  and  normal  stress  bal¬ 
ances  require 

(w,  +  *',)(!  ~hl)  -  2h,(u,  (7) 

P~ahJ{\  +  h]yn 

4-  (2p/(l  +  fi*)]  [u,h2,  -  (uf  +  v,)h,  +  V,}  (8) 

where  n  is  the  surface  tension  coefficient.  At  the  wall,  y  -  0. 

u  ~  K  ,  v  -  0,  (9) 
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represent  the  standard  no  slip  and  no  flux  conditions  Velocities 
at  the  interface  are  related  through  the  kinematic  condition  in  a 
moving  frame. 

v  -  uh,,  y  -  h(z).  (10) 

The  inlet  velocity  profile  is  parabolic,  representing  an  accelera¬ 
tion-free  falling  film,  while  r.  sufficient  and  physically  consistent 
outlet  condition  for  a  solitary  wave  requires  a  zero  slreamwise 
derivative  for  ail  variables  The  variable  replaced  h(z)  as  the 
final  variable  to  be  iteratively  determined  in  the  free-surface 
problem,  and  completes  a  now  well-posed  problem. 

For  each  wave  profile,  a  unique,  nonuniform  finite-difference 
grid  mesh  was  constructed.  The  mesh  for  a  typical  domain  is 
shown  in  Figure  3.  The  particular  wave  shape  determined  the 
grid  spacing  used.  Mesh  refinement  continued  until  no  further 
change  in  either  the  computed  wave  velocity  or  wall  shear  stress 
profile  was  observed.  Of  particular  importance  was  the  concen¬ 
tration  of  cells  near  th  front  and  top  of  the  wave,  since  the 
velocity  fields  change  drastically  in  this  region  due  to  the  large 
interfacial  slope  and  curvature.  For  most  waves,  1,200  cells  of 
dimension  bxby  were  sufficient,  and  produced  grid  Reynolds 
numbers  -  u(z,  y)  hxj v,  ReCy  -  v(z,  y )  oy/v)  of  order  1  in 
the  y  direction,  and  ranging  fr-’m  1  to  100  in  the  streamwise 
direction. 

The  curved  interface  was  accommodated  by  allowing  boun¬ 
dary  cells  to  be  evt  by  the  boundary,  h(z),  thus  reducing  their 
volume.  This  situation  is  illustrated  in  Figure  4.  This  technique 
produced  areas  adjoining  two  boundary  cells,  the  centers  of 
which  were  outside  the  computational  domain.  As  the  stress- 
free  interface  requires  a  zero  normal  derivative  of  the  velocity 
vector  with  respect  to  the  boundary,  h(z),  these  regions  were 
treated  as  inviscid  channels  through  which  all  fluid  leaving  one 
boundary  cell  on  its  shared  side  passed  into  the  neighboring  cell 
through  its  respective  shared  side.  The  total  area  of  these  regions 
represents  less  than  0.1%  of  the  total  domain,  and  had  little 
effect  on  the  results. 

The  equation  set  was  solved  on  a  finite-difference  grid  using  a 
variant  of  the  TEACH-T  code  (Gosman  et  al.,  1969),  incorpo¬ 
rating  the  SIMPLER  pressure-continuity  solution  procedure; 
the  principles  of  this  method  are  described  in  detail  elsewhere 
(Patankar,  1980).  The  domain  includes  regions  of  significant 
streamwise  variation  in  all  variables,  thus  necessitating  an  accu- 
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Figure  3.  Sample  finite-difference  grid. 


Figure  4.  Interface  finite-difference  grid. 

rate  method  of  discretizing  convective  momentum  terms.  The 
simplest  method  of  convective  discretization,  upwind  differenc¬ 
ing  insures  a  reasonably  stable  numerical  solution,  but  intro¬ 
duces  numerical  diffusion  in  regions  of  the  flow  where  stream¬ 
lines  are  oblique  with  respect  to  the  grid  lines  (Raithby,  1976). 
More  importantly,  the  upwide  scheme  lacks  sensitivity  to  cross- 
stream  diffusion  and  source  terms  (Leonard,  1 979),  which  are  of 
tremt  taous  importance  in  the  case  of  a  thin  film.  This  lack  of 
sensitivity  diminishes  the  effects  of  the  y  direction  diffusion  as 
well  as  the  v  velocity  in  the  solution  of  the  streamwise  velocity. 
These  deficiencies  in  the  upwind  and  hybrid  methods  require 
the  use  of  a  QUICK-based  scheme,  which  improves  accuracy  by 
expanding  the  number  of  neighbo.ing  points  included  in  interpo¬ 
lated  values  of  velocity. 

Based  on  Leonard's  (1979)  third-order-accurate  discretiza¬ 
tion  scheme  QUICK,  Pollard  and  Siu  (1982)  developed  the 
QU1CK-ER  (Extended  and  Revised)  method  of  discretizing 
convective  terms.  The  QUICK-ER  method  overcomes  stability 
problems  inherent  in  the  QUICK  procedure  at  the  expense  of 
slower  convergence,  and  is  considered  the  most  satisfactory 
method  of  handling  convec  ive  momentum  terms  (Huang  et  al„ 
1985).  For  application  tc  nonuniform  grids,  a  new  version  of 
QUICK-ER  was  developed.  This  method  follows  the  spirit  of 
the  QUICK-ER  formulation,  but  includes  locally  variable 
weighting  factor  to  account  for  the  nonuniformity  of  the  grid  in 
both  directions.  p  hough  QUICK-ER  schemes  requires  more 
computational  effort  per  iteration  than  upwinding.  particularly 
for  nonuniform  grids,  improvements  in  accuracy  allow  the  use  of 
a  slightly  coarser  grid,  so  total  computational  time  exceeds  that 
required  by  the  upwind  method  by  only  20%. 

The  solution  procedure  began  with  choosing  a  vaiue  for  V, 
and  creating  the  transformed  domain,  given  by  Eq.  2.  The  u 
velocity  field  was  set  to  a  parabolic  profile  everywhere,  and  the  v 
velocity  field  was  set  to  zero.  The  pressure  at  each  z  location  was 
set  to  the  surface  pressure  due  to  curvature.  Updated  velocity 
and  pressure  fields  within  the  wave  were  then  computed  using 
Eqs.  4, 5,  and  6.  Through  interpolation  for  the  velocity  gradients 
in  the  interfacial  shear  stress  balance,  Eq.  7,  streamwise  and 
normal  velocities  in  the  interior  of  the  flow  field  were  used  to 
derive  an  expression  for  the  streamwise  surface  velocity.  Cou¬ 
pled  with  the  kinematic  condition,  Eq.  10,  the  velocities  on  the 
surface  were  known  for  each  iteration.  The  surface  pressure 
computed  from  Eq.  8  was  used  to  determine  me  first  pressure 
value  in  the  interior  of  the  domain  through  the  use  of  parabolic 
interpolation  using  the  surface  pres- .ire  and  two  interior  pres¬ 
sures.  With  the  newly  computed  surface  variables,  the  velocity 
and  pressure  fields  were  updated  until  the  sum  of  residuals  of 
mass  and  momentum  (normalized  by  the  inlet  quantities)  over 
the  domain  was  icss  than  10'  \  This  condition  also  required  the 
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inierfaciai  shear  and  normal  stress  balances  to  be  wit  hi  n  !0  1 1’a 
of  zero,  and  thus,  the  governing  equations  and  boundary  condi¬ 
tions  were  satisfied. 

Upon  convergence  of  the  velocity  and  pressure  fields,  the 
average  pressure  in  the  fiat  outlet  section  was  examined  If  the 
average  pressure  did  not  approach  zero,  as  required  for  a  nonac- 
celcrating  film  surrounding  a  solitary  wave,  a  new  value  of  Vm 
was  chosen  and  the  process  repeated.  The  adjustment  procedure 
for  Vw  was  simple:  If  the  pressure  in  the  outlet  section  was  higher 
than  zero,  the  wave  (wall)  velocity  was  too  high,  since  the  wall 
was  pushing  excess  fluid  through  the  wave,  and  a  positive  pres¬ 
sure  at  the  outlet  was  opposing  this  extra  fluid  in  an  attempt  to 
satisfy  the  mass  balance  for  the  wave. 

Examination  of  the  experimental  data  reveals  that  the  large 
waves  do  not  remain  precisely  constant  in  shape.  Incorporating 
this  unsteady  effect  is  accomplished  through  the  use  of  a  locally 
constant  stretching  parameter,  as  opposed  to  the  globally  con¬ 
stant  value  used  for  the  classical  solitary  wave.  The  domain 
transformation  for  this  case  is  given  by: 

*,  -  *.  +  ~  O  01) 

where  z  is  the  streamwise  coordinate.  In  general,  this  pseudo¬ 
wave  velocity  is 

V„-Vw[\-  *(*,)]  (12) 

where  e(z()  is  an  iteratively  determined  local  stretching  vari¬ 
able,  and  Vw  represents  the  wave  velocity  associated  with  the 
substrate.  The  solitary  wave  case  is  recovered  by  setting 
e(z,)  -  0  for  all  i.  As  before,  we  define  a  new  streamwise  veloc¬ 
ity  component  as 

u(z.y)-u'(x,y)+VJzl)  (13) 

which  a'lows  the  same  governing  equations  and  boundary  condi¬ 
tions,  Eq.  4-10,  to  apply.  The  transformation,  Eq.  12,  introduces 
a  locally  variable  mass  and  momentum  source  due  to  the  evolv¬ 
ing  nature  of  the  wave,  which  does  not  appear  in  a  solitary  wave. 
The  solution  procedure  is  identical  to  that  of  the  purely  solitary 
wave.  The  solution  procedure  is  identical  to  that  of  the  purely 
solitary  wave,  with  the  exception  that  now  a  profile  of  must 

be  specified  instead  of  a  single  value.  When  the  velocity  and 
pressure  fields  have  converged  for  a  given  set  of  V„,  the  wave 
shape  is  adjusted  through  e(z,)  to  meet  two  criteria.  The  base¬ 
line  wave  velocity,  was  adjusted  such  that  the  average  pres¬ 
sure  in  the  fiat  outlet  section  approached  zero,  as  before.  The 
computed  wall  shear  stress  profile  was  then  compared  to  the 
experimental  profile,  and  adjustments  made  to  e(z,)  to  correct 
deviations.  In  this  sense,  the  classical  free-boundary  problem  is 
recovered,  albeit  supplemented  by  experimental  data. 

The  procedure  developed  for  the  solitary  waves  required  an 
average  of  300  iterations  of  the  velocity  and  pressure  fields  to 
converge,  with  an  underrelaxation  factor  of  0.5  used  for  all  vari¬ 
ables.  Between  four  and  eight  adjustments  to  the  solitary  wave 
velocity  were  required  to  produce  a  flow  with  an  average  outlet 
pressure  less  than  IQ"2  Pa.  For  the  quasi-unsteady  case,  the 
same  number  of  iterations  was  required  to  achieve  convergence 
of  the  velocity  and  pressure  fields,  while  the  adjLstment  of  the 
variable  wave  velocity  to  match  wall  shear  stress  data  took  any¬ 
where  from  five  to  20  iterations. 
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Figure  S.  Solitary  wave  streamline  map,  peak /substrate 

«3. 

The  program  was  coded  in  FORTRAN  77,  and  required  2 
MB  of  task  space.  Execution  times  fer  convergence  of  the  veloc¬ 
ity  and  pressure  fields  were  approximately  5  CPU  h  on  a  VAX 
1 1-750  computer. 

Results 

Initial  computations  explored  the  validity  of  neglecting  the 
evolution  of  the  wave  shape.  For  a  pcak/substrate  thickness 
ratio  of  approximately  3,  the  resulting  streamline  map  and  wall 
shear  stress  profile  comparisons  are  shown  in  Figures  5  and  6. 
The  body  force  due  to  gravity,  pgh,  is  presented  to  accommodate 
comparison  with  the  wall  shear  stress  predicted  by  a  parabolic 
velocity  profile.  The  shear  stress  comparison  suggests  that  the 
front  and  back  of  the  wave  arc  accurately  described  as  moving 
undeformed,  but  neglecting  the  evolution  of  the  peak  of  the  wave 
causes  discrepancies  between  computed  and  measured  wall 
shear  stresses.  Figure  7  shows  superposed  traces  of  the  film 
thickness  measured  at  locations  63  mm  apart,  and  clearly  illus¬ 
trates  the  slight  difference  in  speed  between  the  front  and  back 
of  the  wave  as  it  moves  down  the  lube,  suggesting  e(z,)  >  0  in 
this  front  region.  In  order  to  evaluate  the  effect  of  small  changes 
in  velocity  along  the  wave,  the  parameter  e(z,),  given  by  Eq.  12, 
was  varied  by  irial  and  the  effect  on  the  resulting  wall  shear 
stress  comparison  noted.  Figure  8  shows  the  very  slight  degree  of 
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Figure  6.  Solitary  wave  shear  stress  comparison,  peak/ 
substrate  *3. 
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Figure  7.  Evolution  of  nearly  solitary  wave,  peak /sub¬ 
strate  «3. 


variation  in  Vw  along  the  front  of  the  wave,  which  resulted  in 
satisfactory  agreement  between  measured  and  computed  wall 
shear  stress  data  along  the  wave.  To  accomplish  this  match,  the 
wave  velocity  was  increased  gradually  over  the  front  of  the  wave, 
reaching  a  maximum  deviation  of  roughly  10%  directly  under 
the  peak.  The  effect  of  this  variation  on  the  streamlines  and 
shape  of  the  domain  is  small,  but  still  significant,  as  seen  in  com¬ 
paring  Figure  9  with  Figure  6.  While  the  flow  near  the  wall 
appears  quite  sensitive  to  unsteady  effects,  global  flow  patterns 
show  only  small  sensitivity  :o  these  changes.  These  results 
emphasize  the  importance  of  the  interfacial  shape  in  the  compu¬ 
tation  of  the  flow  field. 

Streamline  maps  of  two  additional  waves,  having  peak/sub¬ 
strate  thicknesses  of  approximately  4  and  5,  are  presented  in 
Figures  10  and  1 1.  Streamlines  for  these  larger  waves  were  less 
affected  by  the  transition  from  solitary  to  evolving  waves  than 
the  smaller  wave.  For  these  larger  waves,  the  e(z,)  factors  were 
slightly  larger  than  in  the  previous  wave,  and  the  wall  shear 
stress  comparisons  were  similarly  favorable.  These  streamline 
maps,  in  conjunction  with  Figure  9,  suggest  it  is  reasonable  to 
view  the  waves  as  lumps  of  fluid  overrunning  a  slow  moving  sub¬ 
strate.  Note  that  for  all  three  of  these  waves  of  different  ampli¬ 
tude,  a  well-defined  recirculatory  region  appears  when  viewed  in 
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Figures.  Shear  stress  comparison  (or  evolving  wave, 
peak /substrate  -3. 
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Figure  9.  Streamline  map  for  evolving  wave,  peak  /sub¬ 
strate  -3. 

a  coordinate  system  fixed  on  the  wave.  Of  importance  is  the 
presence  of  large  normal  velocities  near  the  front  and  rear  of  the 
recirculating  region.  Previous  modeling  charactered  all  nor¬ 
mal  velocities  as  small  corrections  to  the  dominant  streamwise 
flow,  a  supposition  now  seen  to  be  inaccurate.  These  normal 
velocities  can  be  expected  to  enhance  transport  by  refreshing  the 
surface  with  fluid  from  the  substrate. 

Surface  velocity  of  the  waves  is  nearly  uniform  over  a  large 
portion  of  the  wave  peak,  while  varying  rapidly  ear  the  front. 
Note  the  presence  of  stagnation  points  in  fron.  of  and  behind 
each  wave  peak.  While  not  physically  important  features  of  the 
flow,  these  points  correspond  to  zero  curvature  of  the  stream- 
wise  velocity,  and  are  unrealizable  for  a  parabolic  velocity  pro¬ 
file. 

Most  previous  modeling  efforts  have  regarded  acceleration 
within  the  wave  as  negligible  compared  to  the  gravitational 
acceleration.  In  the  Kapitza,  or  long-wave,  analysis,  the  inertial 
terms  are  neglected  to  produce  a  linear  hydrodynamic  problem, 
which  yields  the  streamwise  velocity 

u(z.y)  -  (g/»)  [yfi(z)  -  y1/2)  (14) 

For  this  simple  velocity  profile,  it  is  easily  shown  that  the  inertial 
forces,  normalized  with  respect  to  gravity,  are 

( udu/dz  +  vdu/dy)/g  -  (g/y*)y2h(z)/2(dh/dz)  (15) 

which  has  a  maximum  value  at  the  interface,  >>  -  h(z).  given 
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Figure  10.  Streamline  map  tor  evolving  wave,  peak/eub- 
strate  -4. 
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Figure  1 1.  Streamline  map  for  evolving  wave,  peak /sub¬ 
strate  >5. 

by 
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Figure  13.  Convective  acceleration  profile,  peak/sub¬ 
strate  -3. 


gh(z)\dhldz)/U y>)  (16) 

Using  the  physical  properties  given  previously,  this  quantity  (for 
Re  -  880)  becomes  O  (10J)  dh/dz.  For  the  original  lineariza¬ 
tion  of  the  problem  to  be  valid,  Eq.  16  suggests  dh/dz  should  be 
less  than  O  (101).  The  large  waves  discussed  here  have  slopes 
as  great  as  0.07,  in  which  case  inertial  forces  are  predicted  to 
dominate  gravity,  and  the  original  premise  is  violated.  Further 
comparison  of  the  computed  inertial  forces  and  the  long  wave 
predictions,  Eqs.  IS- 16,  must  be  limited  to  the  location  of 
extrema  and  existence  of  inertial  forces  when  the  interfacial 
slope  is  small. 

As  seen  in  the  streamline  maps,  large  accelerations  exist  near 
the  surface,  particularly  near  the  front  and  rear  of  the  recircu¬ 
lating  region.  For  the  wave  having  a  peak/substrate  thickness 
ratio  of  roughly  3  (a  ratio  commonly  seen  in  experimental  data), 
consider  the  three  locations  shown  in  Figure  12.  Near  the  front 
of  the  wave,  maximum  normalized  inertial  forces  (convective 
momentum  terms/gravity),  shown  in  Figure  13,  are  several 
times  greater  than  gravity.  Beneath  the  wave  peak,  roughly  25% 
of  the  flow  is  free  of  acceleration,  a  stagnant  lump  riding  on  the 
substrate.  In  this  region,  the  interfacial  slope  is  nearly  zero,  for 
which  Eq.  16  predicts  no  acceleration  whatsoever.  However,  the 


Figure  12.  Location  of  acceleration  and  velocity  fit  exam¬ 
ples,  peak /substrate  ~3. 


substrate  beneath  the  peak  shows  moderate  deceleration,  with 
the  maximum  occurring  far  from  the  interface,  contradicting 
Eq.  1 5,  which  requires  the  maximum  inertial  force  to  exist  at  the 
free  interface.  Further  toward  the  back  of  the  wave,  the  acceler¬ 
ation  is  again  maximum  at  the  surface,  showing  a  deceleration 
of  the  order  of  gravity,  while  Eq.  16  predicts  a  maximum  accel¬ 
eration  of  the  order  of  10“'  g  for  this  nearly  flat  region.  These 
variations  in  acceleration  in  the  streamwise  direction  support 
the  premise  of  large  waves  pushing  material  from  the  surface  to 
the  substrate,  exchanging  either  heat  or  mass,  and  thereby 
enhancing  transport  properties. 

Acceleration  effects  are  further  illustrated  in  Figure  14;  a 
particle  on  a  surface  streamline  would  undergo  greater  changes 
in  acceleration  than  one  near  the  wall.  The  process  of  transport 
enhancement  is  again  clearly  shown,  as  particles  in  front  of  the 
peak  accelerate  toward  the  mass  of  stagnant  fluid,  exchanging 
heat  or  mass,  and  then  are  forced  to  return  to  the  substrate. 

To  examine  the  suitability  of  various  polynomial  representa¬ 
tions  of  the  streamwise  velocity  profile,  data  from  the  numerical 
experiment  were  compared  to  the  velocity  profile  predicted  by 
Kapitza  analysis, 

u(z,y)  -  2u\z,y  -  h(z )]  \y/h(z)  -  '/2  \y/h(z)Y)  (17) 

where  the  surface  velocity  was  taken  from  the  numerical  experi¬ 
ment.  In  addition,  the  computed  velocity  profile  was  fit  with  a 
least-squares  cubic  polynomial  in  y.  For  illustrative  purposes, 
consider  the  streamwise  locations  within  the  wave  shown  in  Fig¬ 
ure  1 2.  Near  the  front  and  beneath  the  peak,  the  parabolic  fits 
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Figure  14.  Fluid  particle  acceleration,  peek /substrate 
-3. 


AICHE  Journal 


February  1989  Vol.  35,  No.  2 


193 


shown  in  Figures  !  5a  and  1 5b  show  significant  deviations  from 
the  computed  velocities,  while  the  cubic  polynomial  appears  to 
faithfully  represent  the  profile.  Use  of  the  cubic  polynomial 
allows  the  curve  to  fit  the  velocity  gradient  at  the  wall  as  well  as 
the  surface  and  wall  velocity  accommodated  by  the  parabolic  fit 
This  extra  degree  of  freedom  enables  a  nearly  perfect  fit  of  the 


Location  C 


Location  B 


*3. 


data,  and  represents  the  physically  desirable  property  of  match¬ 
ing  wall  shear  stress  as  well  as  surface  velocity.  Near  the  back  of 
the  wave,  where  the  shear  stress  suggests  the  flow  has  small 
acceleration  forces,  the  cubic  and  parabolic  fits  are  shown  in 
Figure  I  Sc  to  coincide  more  closely  than  at  other  locations 

In  the  cases  of  both  solitary  and  evolving  waves,  the  wave 
velocity  was  computed  by  choosing  K.  such  that  the  average 
outlet  pressure  was  close  to  zero.  The  dependence  of  average 
outlet  pressure  on  the  wave  velocity  is  shown  in  Figure  16  for  the 
solitary  wave  with  a  pcak/substrate  ratio  of  approximately 
three.  This  nearly  linear  dependence  allowed  quick  convergence 
to  yw.  For  each  wave,  the  substrate  or  base  wave  velocity  Vw  was 
compared  to  that  determined  by  dividing  the  length  between  the 
upper  and  lower  film  thickness  probes  by  the  passage  time  of  the 
wave  peak.  In  cases  where  the  wave  traveled  between  the  probes 
with  only  slight  deformation,  this  comparison  showed  the  com¬ 
puted  values  to  be  within  roughiy  10%  of  the  experimental  ones. 
For  those  waves  evolving  rapidly,  this  crude  comparison  pro¬ 
duced  less  favorable  results,  with  errors  as  high  as  30%  for  the 
wave  and  peak/substrate  of  approximately  5. 

Figure  17  shows  the  values  of  non-dimensional  wave  velocity 
determined  by  simulating  flow  in  various  solitary  waves  at  dif¬ 
ferent  flow  rates.  The  hydrodynamic  character  of  these  waves 
was  similar  to  those  presented  in  detail:  of  primary  interest  was 
the  lack  of  correlation  of  wave  velocity  to  peak/substrate  thick¬ 
ness,  in  agreement  with  the  experimental  findings  of  Zabaras 
(1985). 

Conclusions 

The  interface  of  a  falling  liquid  film  consists  of  a  random 
array  of  waves  of  varying  amplitude,  length,  and  velocity,  some 
isolated  and  some  overlapping.  Even  at  moderate  Reynolds 
numbers,  these  waves  display  amplitudes  that  are  two  to  five 
times  the  substrate  thickness.  Three  typical  isolated  waves 
obtained  from  measurements  of  the  time  traces  of  the  film  thick¬ 
ness  were  selected  as  computational  domains.  A  method  was 
developed  to  solve  the  Navier-Stokes  equations  for  this  free-sur- 
face  problem  which  yields  both  velocity  and  pressure  fields  in 
the  wave  and  wave  velocity.  The  measured  wall  shear  stresses 
and  wave  velocities  were  in  reasonable  agreement  with  those 
determined  from  the  numerical  solutions.  The  results  of  these 
computations  confirm  certain  earlier  speculations  on  the  me- 


Flgure  16.  Average  outlet  pressure  dependence  of  wave 
velocity,  peak/substrote  *  3. 
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Figure  17.  Numerical  prediction*  of  wave  velocity. 

chanics  of  wavy  film  flow  and  point  to  the  inadequacies  of  oth¬ 
ers. 

The  large  wave  moves  rapidly  over  a  slow  substrate,  the  bulk 
of  the  liquid  being  carried  in  the  wave;  this  liquid  is  nearly  sta¬ 
tionary  in  a  coordinate  system  moving  with  the  wave.  The  inter¬ 
action  between  the  wave  and  substrate  causes  the  acceleration  of 
fluid  at  the  front  of  the  wave  from  the  substrate  in  the  peak.  The 
fluid  then  decelerates  as  it  passes  out  the  wave  back,  generating 
a  close  recirculation  region  in  the  wave  whose  size  depends  on 
the  wave  amplitude.  This  process  of  recirculation  would  appear 
to  account  for  the  enhanced  rates  of  heat  and  mass  transfer 
known  to  exist  in  wavy  films.  This  picture  of  the  waves  is  in 
rough  accord  with  the  speculations  of  Dukler  (1977). 

The  presence  of  secondary  flows  can  be  expected  to  enhance 
rates  of  interfacial  mass  transfer.  Lcvich  (1962)  showed  that 
velocities  normal  to  wall  as  small  as  1%  of  the  streamwise  veloc¬ 
ity  produces  a  1 5%  increase  in  the  rate  of  mass  transfer  to  the 
interface.  As  the  present  work  shows  that  normal  velocities  may 
be  ten  times  this  large,  we  expect  substantial  enhancement. 

Most  previous  studies  of  wavy  film  flows,  starting  with  the 
classical  work  of  Kapitza  (1964),  are  based  on  the  parabolic 
streamwise  velocity  profile,  which  is  presumed  to  exist  at  all 
positions  along  the  wave.  This  numerical  experiment  points  to 
the  inadequacy  of  that  assumption  for  these  large  waves,  and 
shows  that  a  cubic  profile  faithfully  reproduces  the  streamwise 
velocity  gradient  at  all  axial  locations.  Attempts  to  develop  evo¬ 
lutionary  equations  for  this  flow  should  use  the  higher  order 
velocity  profile  to  insure  that  the  accelerations  which  exist  can 
be  incorporated  into  the  model. 

Most  theoretical  models  have  shown  that  the  wave  velocity 
and  amplitude  are  uniquely  related.  Careful  measurements  by 
Zabaras  (1985)  show  that  the  wave  amplitude  varies  widely 
with  both  amplitude  and  wave  length.  These  numerical  studies 
show  that  for  a  given  amplitude,  the  computed  wave  velocity  is 
sensitive  to  relatively  small  changes  in  wave  shape,  rationalizing 
the  experimental  observations  and  raising  questions  as  to  the 
usefulness  of  some  of  the  simpler  models. 

The  use  of  a  single  wavelength  as  a  scaling  parameter  for  the 
streamwise  hydrodynamics  is  shown  to  be  incorrect.  The  present 
computations  show  that  strong  differences  in  accelerations  exist 
at  various  positions  along  the  wave.  Each  of  these  domains 
should  be  analyzed  separately  if  integral  equations  are  to  be 
used,  as  suggested  by  Maron  et  al.  (1985).  The  use  of  evolution¬ 
ary  equations  can  eliminate  this  problem. 
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Notation 

D  -  liquid  -liquid  diffusion  coefficient,  m'/t 
g  -  acceleration  of  gravity,  m/s’ 

P  -  pressure,  N/mJ 

Q  -  liquid  film  flow  rale  per  unit  perimeter,  m’/s 
Re  -  film  Reynolds  number,  4 Q/e 
Sc  -  Schmidt  number,  e/D 
u  -  local  streamwise  velocity,  m/s 
v  -  local  velocity  normal  to  boundary,  m/s 
V.  -  solitary  wave  velocity,  m/s 
x  -  axial  coordinate  in  lab  frame,  m 
y  -  coordinate  normal  to  boundary,  m 
i  -  coordinate  fixed  on  wave,  m 
P  -  liquid  density,  kg/m1 
v  -  liquid  kinematic  viscosity,  m’/s 
A  -  Laplacian  operator,  dfdx1  +  if/d/ 
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Phase  Plane  and  Bifurcation  Analysis  of  Thin 
Wavy  Films  under  Shear 


A  long-wave  equation  for  film  thickness  as  a  function  of  position  is 
derived  for  a  general  case  incorporating  viscous,  surface  tension,  and 
interfacial  shear  effects.  The  derivation  considers  both  the  parabolic 
and  the  power-law  velocity  profiles.  The  analysis  is  aimed  at  revealing 
the  wave  velocity  that  induces  infinitely  long  (homoclinic)  periods  as 
well  as  substrate  thickness  and  wave  peak  amplitude.  Phase  plane 
analysis  shows  that  at  Re  »  1,  due  to  time-scale  separation,  the  homo¬ 
clinic  velocity  is  near  that  at  the  Hopf  bifurcation.  That  enables  analyti¬ 
cal  derivation  of  the  wave  characteristics. 

Comparison  with  experimental  results  in  the  range  of  ffe-3t0-3,100 
with  countercurrent  gas  flow,  shows  encouraging  agreement.  At  very 
high  Re  the  wave  velocity  suggests  the  onset  of  turbulence,  in  agree¬ 
ment  with  theory.  Phase  plane  analysis  predicts  also  that  the  wave 
shape  consists  of  a  simple  peak  with  a  steep  front,  with  short  waves 
riding  on  the  main  wave  at  low  Re. 


M.  Sheintuch,  A.  E.  Dukler 

Chemical  Engineering  Department 
University  ot  Houston 
Houston.  TX  77204 


Introduction 

The  prediction  of  wave  characteristics  of  falling  liquid  films 
has  been  the  subject  of  numerous  investigations  (Dukler,  1977). 
Most  studies  have  focused  on  free-falling  films  at  low  Reynolds 
number  (Kapitza  and  Kapitza,  1949;  Alekseenko  et  al.,  1985) 
and  have  used  linear  stability  analysis  to  determine  the  range  of 
unstable  wave  velocities  and  the  mode  of  the  fastest  growing 
wave.  The  film  is  assumed  to  be  sinusoidal,  oscillating  about  its 
mean  value  or  even  somewhat  perturbed  from  a  sinusoidal 
shape.  This  type  of  analysis  predicts  the  wavelength  and  velocity 
at  the  conditions  of  wave  inception.  The  long-wave  approxima¬ 
tion  (Benney,  1966)  considers  the  limit  of  very  long  and  shallow 
waves  that  reach  constant  shape  in  the  frame  of  a  moving  coor¬ 
dinate.  Again,  the  wave  evolves  from  the  mean  film  thickness. 
Recent  nonlinear  analysis  searched  for  wave  velocities  that 
induce  infinitely  long  periodic  solutions  (homoclinic  oTbits). 
These  studies  were  limited  to  falling  films  at  low  Reynolds  num¬ 
bers  (Pumir  et  al.,  1983)  or  used  a  questionable  approximation 
for  the  flow  condition  (Needham  and  Merkin,  1984).  Another 
approach,  applied  at  high  Re,  is  to  approximate  the  wave  shape 
by  assuming  a  sequence  of  characteristic  velocity  profiles  at  dif- 
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ferent  locations  along  the  wave  (Maron  et  al.,  1985).  Finally, 
numerical  solutions  of  the  evolution  equation  (Bach  and  Villad- 
sen,  1984)  did  show  the  existence  of  constant -shape  solitary 
waves.  Even  at  Reynolds  number  of  order  10  such  computations 
are  highly  time  consuming,  and  problems  of  numerical  instabil¬ 
ity  are  encountered. 

The  purpose  of  this  analysts  is  to  present  asymptotic  solutions 
for  wave  velocities  as  well  as  certain  dimensional  characteristics 
of  falling  liquid  films  experiencing  interfacial  shear  induced  by 
gas  flow.  For  this  purpose  methods  of  nonlinear  analysis  are 
used.  Experiments  have  shown  that  at  positions  well  below  the 
plane  of  wave  inception,  two  classes  of  waves  exist  on  the  surface 
(Dukler,  1977).  Large  waves  can  be  observed  riding  on  a  thin 
substrate,  with  the  ratio  of  peak  wave  height  to  substrate  thick¬ 
ness  ranging  from  2  to  4.  A  second  class  of  waves  that  are  capil¬ 
lary  in  nature  and  of  much  smaller  amplitude  than  the  first  is 
also  present.  These  ride  on  the  substrate  and  are  sometimes  seen 
on  the  trailing  edge  of  the  large  waves.  In  the  presence  of  gas 
flow  parallel  to  the  mean  surface,  interfacial  shear  is  generated 
and  the  wave  structure  changes  but  the  two  wave  class  structure 
can  still  be  seen.  At  sufficiently  high  countercurrent  gas  flow 
rates  a  condition  of  flooding  can  be  reached  where  part  of  the 
liquid  introduced  on  the  vertical  surface  flows  upward. 

For  complex  problems  such  as  this  one,  two  rather  divergent 
approaches  are  available.  At  low  Reynolds  -•■'mbers  an  analyti- 
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cal  method  is  possible  in  the  vicinity  of  certain  singular  points 
that  can  be  shown  to  exist  (Chang,  1986).  For  high  Re  a  dif¬ 
ferent  approach  is  possible  if  it  can  be  shown  that  the  system  is 
characterized  by  two  or  more  widely  different  time  scales.  This 
paper  takes  the  Utter  avenue,  which  has  been  largely  unex¬ 
plored.  First  we  derive  the  integral  momentum  equation  for  the 
film  under  conditions  of  constant  interfacial  shear  along  the 
wave.  Surface  tension  effects  are  incorporated  and  we  assume 
alternatively  a  parabolic  (laminar)  or  a  power  law  (turbulent) 
profile.  A  complete  analysis  of  the  parabolic  profile  case  is 
developed  including  interfacial  shear  while  for  turbulent  flow 
the  solution  is  explored  only  for  free-falling  films  in  the  absence 
of  gas  flow. 

The  analysis  strategy  is  then  outlined  and  this  is  followed  by 
the  development  of  the  solution  to  the  flow  equations  for  the  two 
profiles,  neglecting  surface  tension.  For  the  parabolic  profile  it  is 
shown  that  the  shear  makes  no  qualitative  changes  in  the  behav¬ 
ior  of  the  system,  having  only  a  quantitative  effect  on  the  inter- 
facial  characteristics.  As  long  as  flooding  is  not  induced  the 
results  are  similar  to  those  of  the  free-falling  case.  This  same 
conclusion  was  reached  by  Zabaras  and  Dukler  (1988)  based  on 
experimental  observations.  Turbulent  flow  may  significantly  re¬ 
duce  the  wave  velocity  to  a  value  approaching  the  mean  liquid 
film  velocity.  Then  the  role  of  surface  tension  is  explored.  In  the 
absence  of  surface  tension  an  ordinary  differential  equation  of 
the  second  order  is  generated  and  it  is  shown  that  the  smooth 
film  (Nusselt)  solution  cannot  be  a  saddle  point.  When  surface 
tension  is  added,  a  third-order  ordinary  differential  equation  is 
generated  and  every  solution  may  be  the  saddle  point  and  thus 
can  be  the  source  of  the  homoclinic  orbit.  The  analysis  of  this 
complete  form  of  the  problem  reveals  that  capillary  waves  can 
exist  along  the  wave  and  substrate  where  the  profile  is  flat,  in 
addition  to  the  large,  long  waves  characterized  by  the  homo¬ 
clinic  orbit;  this  too  is  in  accord  with  experiment. 

Finally,  we  compare  the  results  with  experiments  of  Zabaras 
and  Dukler  (1988). 


With  integration  in  they  direction,  Eq  I  yields 
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with  the  local  flow  rate 
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Assume  a  quadratic  velocity  profile  that  satisfies  the  conditions. 
u  -  0  at  y  -  0;  t,/m  -  dujdy  at  y  -  h;  and  Q{h)  -  u  dy.  This 
yields  the  following  distribution  function 
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with  t)  -  y/h.  Clearly,  this  equation  also  describes  the  distribu¬ 
tion  in  the  undisturbed  film  (steady  state)  as  well.  With  this  dis¬ 
tribution  the  following  integrals  can  be  evaluated. 
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From  Eq.  6  the  second  x  derivative  of  Q  can  be  evaluated. 
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For  the  parabolic  profile  the  last  integral  can  then  be  evalu¬ 
ated. 


The  Film  Equation 

The  Navier-Stokes  equations  are  integrated  in  the  direction  y 
perpendicular  to  the  wall,  using  both  parabolic  and  power-law 
velocity  profiles.  The  work  may  be  extended  to  higher  polyno¬ 
mial  velocity  profiles  by  numerical  methods.  We  then  apply  a 
moving  coordinate  frame  and  arrive  at  one  third-order  equation 
for  the  film  thickness. 
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Substituting  these  quantities  into  Eq.  5  gives  the  following  equa¬ 
tion  for  the  case  of  constant  shear  independent  of  x. 


Parabolic  laminar  profile 

The  film  flowing  downward  under  interfacial  shear,  t„  is 
described  by 
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with  u^dy  defined  in  Eq.  10.  The  integral  of  the  continuity 
equation  is 


subject  to 


u-o-0  at  y-0  (3) 

u,  -  r ,/m,  v  -  h,  +  uh„  P  -  P0  -  oh „  at  y  -  h  (4) 
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Thus  Eqs.  1 1  and  1 2  are  two  equations  in  h  and  Q  as  dynamic 
variables  and  x  and  1  as  independent  variables.  In  a  coordinate 
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frame  moving  with  the  wave  velocity,  c,  set 


to  be  independent  of  the  position  along  the  wave  I  and  the  coeffi¬ 
cients  arc: 


tf-f,  £-x-rr  (H) 

and  the  mass  conservation  condition  then  becomes 

h,-cht  +  Q(-0  (14) 

The  momentum  balance,  Eq.  1 1 ,  can  be  modified  by  substitut¬ 
ing  Q,  —  cQ(  for  Q„  and  for  d’/dx".  We  make  now  two 

important  assumptions: 

1.  The  waves  propagate  without  change  in  shape;  thus  in  a 
moving  coordinate  system,  Q,  -  h,  -  0.  Integration  of  Eq.  14 
yields 
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Q-ch-K  (15) 

where  Q  is  the  volumetric  flow  rate  at  any  location  £  along  the 
wave. 

2.  The  integration  constant,  K  should  satisfy  the  steady  state 
solution  with  h  -  ha  and  Q  -  Qr,  the  liquid  feed  rate.  Thus  K  - 
Qf  —  ch0.  By  integrating  the  velocity  distribution  over  the 
steady  state  film  thickness  we  find  that 

Qf-~d+T)  (16) 

with  the  dimensionless  shear  T  -  3r,/2pgA„.  Note  that  Eq.  16 
may  have  more  than  one  real  solution  for  h„  given  Qr  and  r„  as 
discussed  in  detail  by  Maron  and  Dukler  (1984).  Substituting  K 
calculated  from  the  steady  state  solution,  Eq.  15  becomes 

6-^0  +  r)[l  +n(H-  1)],  H-~,  n-^2  (17) 
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where  n  is  the  wave  velocity  made  dimensionless  in  respect  to  the 
average  velocity.  Substitution  of  Eq.  17  into  the  momentum  bal¬ 
ance,  Eq.  1 1 ,  yields  an  ordinary  third-order  differential  equation 
describing  the  change  of  dimensionless  film  thickness,  H,  with 
dimensionless  length  /  -  £/£.  The  length  scale,  L,  and  dimen¬ 
sionless  parameters  are  defined  as 
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The  resulting  equation  is 
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where 
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In  this  development  the  dimensionless  shear  stress  r,  is  assumed 


These  equations  apply  for  downflow  of  liquid  with  counter- 
current  gas  flow  (g  >  0,  t(  <  0,  - 1  s  T  <  0)  or  cocurrent  flow 
( g  >  0,  r,  >  0,  T  >  0)  as  well  as  for  upflow  where  x  or  /  is  positive 
in  the  flow  direction  (g  <  0,  -t,  >  0,  T  <  -1);  the  inequality, 
T  s  - 1,  follows  from  the  requirement  that  Qf  be  positive  in  Eq. 
16. 

A  schematic  representation  of  the  relation  between  T  and  ha 
as  deduced  from  Eq.  16  appears  in  Figure  1  for  constant  QF- 
Downflow  with  countercurrent  shear  takes  place  for  —  1  £  T  <  0 
(band  B  in  the  figure).  However  in  the  range  —  3/«  <  T  <0  the 
velocities  are  directed  uniformly  downward  at  all  locations  in 
the  film.  For  T  -  -J/«  the  profile  is  symmetric,  with  zero  velocity 
at  the  interface,  and  when  -  1  <  T  <  —  ’/»,  downward-directed 
velocities  exist  near  the  wall,  with  upward  velocities  near  the 
interface.  Band  A  in  Figure  1  represents  concurrent  downward 
flow  while  band  C  displays  the  behavior  for  concurrent  upward 
flow.  In  these  bands,  for  the  range  - 1 .5  <  T  <  -  1 ,  the  velocity 
profile  is  again  not  monotonic.  For  more  negative  values  of  T 
uniformly  upward  flow  exists  in  the  film.  At  T  -  —1.5  the  wall 
shear  is  zero  and  the  upward  motion  of  the  liquid  is  driven  by  the 
interfacial  shear  only.  The  existence  of  these  different  patterns 
of  velocity  distribution  was  discussed  by  Maron  and  Dukler 
(1984)  and  their  existence  confirmed  experimentally  by  Zaba- 
ras  et  al.  (1986). 
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Figure  1.  Solutions  to  steady  state  film  thickness  equa¬ 
tion. 
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Turbulent  profile  in  the  absence  of  shear 

For  turbulent  (low.  the  film  velocity  profile  can  be  expressed 
in  a  power  law  form: 
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We  find 
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The  integral  momentum  equation  then  becomes: 
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For  turbulent  flow  the  wall  shear  stress  t„  cannot  be  approxi¬ 
mated  from  the  velocity  profile  and  correlations  are  needed.  As 
shown  below,  this  has  no  influence  on  the  dimensionless  wave 
velocity  but  does  change  the  mean  film  thickness  and  conse¬ 
quently  the  dimensional  velocity.  From  the  Blasius  relation 
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where  the  friction  factor  at  the  wall,  X,  is  estimated  from  single¬ 
phase  flow  models. 
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When  this  relation  is  substituted  into  Eq.  24  for  the  flat  film,  we 
find 
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In  the  application  of  Eq.  23  at  high  flow  rates  the  value  of  a 
can  be  expected  to  vary  along  the  wave.  At  the  substrate  the 
local  flow  rate  is  low,  turbulence  is  suppressed,  and  the  qua¬ 
dratic  law  will  apply.  In  the  thicker  portions  of  the  wave  near  the 
crest,  turbulent  behavior  will  result  in  increased  values  of  a. 
Methods  for  accounting  for  this  flow  direction  variation  in  a 
have  not  yet  been  fully  implemented.  As  a  first  approximation 
we  use  a  constant  value  of  a  for  the  dynamic  analysis  of  the 
momentum  equation. 

Strategy  of  Dynamic  Analysis 

The  dynamic  equation  for  the  film  thickness,  Eq.  19,  contains 
four  parameters.  Re,  W,  n,  and  T,  which  are  to  be  specified  as 
input  variables.  Even  though  the  slopes  of  the  waves  are  known 
to  be  small  it  is  not  possible  to  simplify  the  equation  by  neglect¬ 
ing  the  Hu  and  Hm  terms  since  these  quantities  are  necessary  to 
the  stability  and  dynamic  analysis  that  follows.  Define  the  fol¬ 
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When  0  is  large  the  inertial  forces  are  large  compared  to  the  sur¬ 
face  forces  and  the  first  term  in  Eq.  19  can  be  neglected. 

The  steady  state  solution  of  Eq.  19  is  the  space-independent 
(H,  -  H„  -  H„i  -  0)  solution  of  C0(H)  -  0.  The  qualitative 
nature  of  these  solutions  is  shown  in  Figure  2  for  different 
ranges  of  T.  The  unperturbed  film  (H  -  I)  is  always  a  solution 
but  there  is  always  a  range  of  velocities  n  for  which  it  is  not  sta¬ 
ble.  The  other  two  solutions,  for  any  T,  form  a  parabolic  branch 
n  -  H1/(  1  +  T)  +  (1  +  H)  which  intersects  H  -  1  at  n  -  (3  + 
27*)/(l  +  T)  and  acquires  a  turning  point  at  H  -  —(1  +  7')/2, 
/i  —  ( 3  -  7')/4.  The  parabolas  for  upflow.  Figure  2a,  and  down¬ 
flow,  Figure  2b,  open  in  opposite  directions.  The  range  of  values 
of  n  over  which  the  H  -  1  solution  is  unstable  can  be  determined 
from  linear  stability  analysis.  Specifically,  there  will  usually 
exist  a  (Hopf)  bifurcation  point,  which  is  the  value  of  n  at  which 
a  transition  to  periodic  behavior  takes  place  and  oscillations  are 
observed  in  the  frame  of  a  moving  coordinate;  that  is,  waves 
propagate  along  the  film  in  the  physical  system.  The  H-n  dia¬ 
gram  for  T  -  0  showing  the  nature  of  the  various  states  as  deter¬ 
mined  by  linear  stability  analysis  is  shown  in  Figure  3.  As  the 
waves  grow,  their  velocities  vary  in  the  direction  of  lower  n 
(more  unstable).  The  waves  form  a  family  of  closed  curves  (limit 
cycles)  in  the  phase  plane  (H,  vs.  H)  or  phase  space  (//„,  H,,  H) 
and  these  closed  curves  increase  in  area  as  the  waves  grow  in 
amplitude  along  the  film  in  the  coordinate  direction.  These 
waves  have  a  peak  such  that  H  >  1  and  a  substrate  thickness 
H  <  1 .  As  the  wave  grows  the  substrate  becomes  smaller.  At  the 
condition  where  the  substrate  thickness  is  equal  to  that  of  the 
saddle  point  the  wave  can  grow  no  more  and  its  velocity  will  not 
change.  In  the  phase  plane  the  largest  possible  wave  will  be 
formed  when  the  limit  cycle  hits  a  saddle  point  forming  a  homo¬ 
clinic  curve.  This  curve  provides  the  asymptotic  amplitude  of  the 
wave  and  the  corresponding  velocity,  n,.  Further  change  in  n 
cannot  take  place  since  this  results  in  the  elimination  of  the 
oscillations.  The  wavelength  of  this  form  is  infinite  since  the 
motion  takes  place  in  the  vicinity  of  the  saddle  point,  which  is  a 
steady  state,  and  thus  is  infinitely  long.  The  saddle  point  satis¬ 
fies  the  steady  state  solution  C0(H})  -  0. 

The  purpose  of  the  nonlinear  analysis  is  to  determine  the  loca¬ 
tion,  n„  and  form  of  the  homoclinic  curve  (also  termed  saddle- 
loop  bifurcation).  That  is  possible  analytically  only  in  two 
cases: 

1.  In  the  small  neighborhood  of  higher  order  singularities 
such  as  when  Hopf  and  saddle-node  bifurcations  coalesce.  Then 


Figure  2.  Variation  of  steady  state  film  thickness  with  n. 

a.  Downflow;  b,  c.  Upflow 

Coordinates  for  limit  points  in  (b)  and  (c)  same  as  in  <a) 
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Figure  3.  Changes  in  the  nature  of  the  steady  states. 

Determined  from  linear  itability  analysis  for  T  -  0  and  negligible 
surface  forces 

a.  n  -  nH;  b  r<nH 


2.  Exchange  of  stability  from  H  -  I  to  the  H  *  1  solution 
occurs  when 

|y|_  _  C„„(l)-^(l)  (34) 

tjC2  <jCj  dH 

yielding  n  -  3  when  T  -  0,  that  is,  at  the  intersection  point  of  the 
two  branches.  These  transitions  are  observed  in  Figure  3  where 
the  unstable  branches  are  marked  by  broken  lines. 

Oscillations  (waves)  exist  for  a  certain  range  of  n  as  shown  in 
Eq.  33.  We  show  now  that  when  t2  -  27  (I  +  T)*/ Re1  — -  0  the 
wave  velocity  at  the  homoclinic  orbit  is  very  close  to  that  at  the 
Hopf  {nt  —i >  nH).  Note  that  <2  does  not  affect  the  stability  bound¬ 
aries,  Eq.  32,  but  it  affects  the  behavior  in  the  phase  plane.  From 
Eq.  30,  with  T  -  0  and  when  the  small  ^CjW2  term  is  ignored, 
the  direction  of  the  motion  (slope  in  the  H  -  H,  plane)  is 


the  oscillations  are  infinitesimally  small  and  can  be  described  by 
perturbation  of  the  linear  solution. 

2.  When  the  system  is  characterized  by  different  time 
scales. 

We  capitalize  here  on  the  latter  property. 

The  dynamic  analysis  of  Eq.  19  requires  working  in  the  three- 
dimensional  phase  space,  H,  H,,  H„.  However  many  of  the  es¬ 
sential  dynamic  features  of  the  system  can  be  discerned  by 
analyzing  the  phase  plane  H,  H,.  Equation  19  reduces  to  this 
problem  for  high  Reynolds  numbers  when  0  — *  <».  Thus  the  pro¬ 
cedure  will  be  to  first  analyze  the  simplified  system  in  the  phase 
plane  for  both  laminar  and  turbulent  flow.  Then  the  analysis 
will  be  generalized  for  cases  where  fi  is  not  large. 

Laminar  Film  Flow  for  /?  — ►  °° 

In  the  case  W  -  0,  Eq.  19  can  be  written  in  the  form 


-  {H  -  H4)(H  +  Ht)w  -  (//  -  [)(H  -  H1)(H  -  H ,) 
dw  5 


dH 


wH-i-H+H,) 


f{H. ») 
<2wCAH) 


(35) 


where  Hu  are  the  roots  of  C0(H)  -  0  with  //,  and  //,  defined 
from  Eq.  20  after  setting  T  -  0. 

//  -  -0.5  ±  Vn  -  0.75 


H4-  & 


n  -  1 
n 

3(n  -  1) 
n 


(36) 


//,  -  w  (29) 

t2C2w,  -  -(1  +  T)2C,w  -  C0  +  <iC,w>  -  f(H,  w)  (30) 

Linear  stability  analysis  requires  the  determination  of  the  eigen¬ 
values  of  the  Jacobian  matrix  d(H',//«2C2)/d(//,  w)  at  the 
steady  state  H  -  1,  w  -  0.  The  matrix  is  then 


Since  e2  is  small,  dwjdH  — ►  »  everywhere  in  the  phase  plane 
except  when  the  numerator  vanishes.  At  this  condition  of/  -  0 

{37) 

j  (ff  -  H4 )(H  +  H,) 


J  - 


1 


h_  A 

\«2Cj 


(31) 


and  the  state  may  destabilize  by  one  of  two  ways: 

1 .  A  Hopf  bifurcation  to  periodic  behavior  occurs  at 


trJ- 


f2C2 


U  +  T)2C,(1) 
<jC2(D 


-0  |/|>0 


(32) 


The  state  is  unstable  for  C,(l)/C2(l)  <  0,  or  from  Eq.  20  with 

r-o. 


3  6  S 


(33) 


The  shape  of  the  w  vs.  H  curves  resulting  from  this  expression 
are  dependent  on  the  velocity,  rt,  since  H,_s  are  functions  of  n. 
Figure  4a  shows  the  solution  of  Eq.  37  at  n  -  nH  -  1 .689,  Eq.  33. 
H  -  1  is  one  branch  of  the  solution  for  all  w  since  Ht  -  1  at  n  - 
nH.  The  second  branch  crosses  this  H  -  1  branch  and  intersects 
the  abscissa  at  H}.  We  have  already  shown  that  for  any  point  in 
the  phase  plane  the  trajectory  must  be  vertical  except  near  the 
/  -  0  curves.  Thus  all  trajectories  starting  at  H  <  1  flow  toward 
the  second  branch  except  for  those  that  originate  in  the  immedi¬ 
ate  vicinity  of  the  H  -  1  line.  Along  the  curves  the  trajectory 
must  be  horizontal.  As  the/-  0  lines  are  approached  the  trajec¬ 
tory  must  move  along  the/  -  0  curve  and  increasingly  close  to  it 
in  direction.  This  direction  is  indicated  in  the  usual  way  by 
heavy  arrows  on  the  solution  branches.  The  sign  of  w,  is  an  indi¬ 
cator  of  the  direction  of  the  vertical  trajectories.  As  seen  in  Eq. 
30,  this  direction  changes  as /(//,  w)  passes  through  zero.  Simi¬ 
larly,  the  sign  of  H,  is  an  indicator  of  the  direction  of  horizontal 
motion.  Eq,  29  shows  that  this  direction  changes  as  w  passes 
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Figure  4.  Phase-plane  analysis  of  two-variable  model. 

Heavy  lines./ -  0:  thin  lines,  direction  of  motion 


a.  n  <  n„:  b.  n„  <  n  <  n, 
c.  n  -  n,  for  Re  —  «•;  d.  n  -  n, 

ihrough  zero.  Equation  35  shows  that  trajectories  cannot  cross 
the  H  -  H}  boundary  since  at  this  point  Cj(H)  -  0  and  a  discon¬ 
tinuous  change  in  the  direction  of  motion  occurs.  When  this 
oegins  to  happen  the  surface  tension  terms  neglected  here 
oecome  important.  Thus  waves  that  initiate  from  the  smooth 
film  cannot  achieve  an  amplitude  greater  than  Hy 
Solutions  of  Eq.  37  for  n  <  n„and  n  >  n„  appear  in  Figures  4b 
ind  4d.  In  the  latter  case  the  steady  state,  H  -  1 ,  is  stable  and  all 
trajectories  in  its  vicinity  are  attracted  there.  For  n  <  nH  in  Fig¬ 
ure  4b,  the  steady  state  at  H  -  1  is  unstable  and  ail  trajectories 
:scape  to  infinity.  Of  particular  interest  is  the  case  where  n  — 
t„,  illustrated  in  Figure  4c.  The  Hopf  theorem  indicates  that  a 
limit  cycle  surrounds  the  unstable  state,  H  -  1,  near  the  condi- 
Mon  for  a  Hopf  bifurcation.  Thus  the  situation  appears  as  shown 
n  Figure  4c  for  n  somewhat  less  than  n„.  Periodic  oscillations 
wind  out  in  a  spiral  from  its  origin  at  H  -  1  and  gradually  grow 
in  amplitude  until  they  disappear  in  a  homoclinic  orbit  at  Hy 
Now  it  is  possible  to  deduce  that  the  value  of  n,  at  which  this 
akes  place  must  be  negligibly  different  from  nH  when  «j  is  small. 
Under  these  conditions  the  trajectories  are  almost  vertical,  and 
the  limit  cycle  must  be  narrow  in  the  //direction  and  long  in  the 
v  direction.  If  the  two  branches  of  the/  -  0  curve  are  far  apart, 
ts  they  are  when  n  is  significantly  lower  than  nH,  the  trajectories 
can  readily  escape  to  infinity.  As  <2  increases,  the  trajectories 
continue  to  flow  toward  the  left  branch  of  the/  -  0  curve  but 
tow  they  can  travel  a  short  distance,  1 2,  away  from  these  curves. 
The  separation  distance  between  the  two /-  0  curves  increases 
as  n,  —  nH.  Thus  for  the  trajectory  to  move  from  H  -  1  on  one 
/ranch  of/-  0  to  //,  on  the  other  requires  that  n,  -  nH  be  of  the 
irder  of  €},  a  small  number  under  the  experimental  conditions  of 
interest  here.  Thus  we  conclude  that  n,  -  nH. 

Pictures  of  the  developing  waves  can  be  arrived  at  by  ana- 
yzing  the  trajectory  in  the  phase  plane,  as  shown  in  Figure  5. 
The  initially  flat  interface  (a)  displays  small  sinusoidal  waves  at 
the  point  of  wave  inception  (b)  with  the  wave  velocity,  n  -  nH  - 
1 .689,  Eq.  33,  at  this  condition.  But  the  condition  is  unstable  and 
he  spiral  unwinds  around  H  -  l  to  become  a  homoclinic  orbit  at 
i  -  n,  in  close  proximity  to  nM.  In  its  homoclinic  orbit  the  mini¬ 
mum  or  substrate  wave  amplitude  is  //,  -  0  47  (Eq.  36  for  n  - 
iH).  H  then  increases  rapidly  at  the  front  of  the  wave  with  high 


Figure  S.  Qualitative  sketch  ot  wave  profile. 


slope,  w.  As  it  approaches  its  peak  amplitude  the  slope  drops  to 
ze;o.  The  peak  cannot  exceed  //,  -  1.25.  Then  the  sign  of  the 
slope  changes  and  along  the  back  the  wave  tapers  more  slowly 
down  to  the  substrate  thickness,  Hy 
The  existence  of  countercurrent  shear  does  not  change  this 
qualitative  picture,  modifying  only  the  wave  characteristics. 
Steady  states  are  linearly  unstable  when  C,(l)/Cj(l)  <  Oor 


/3  6  s  S  I  22s  41? 

(2  +  2]<"<"H-5  +  25  +  TV1+lT+'9?  (3> 

where  s-T/(l  +  T).  The  upper  root  of  C,(l)  -  0  is  n„and  this 
is  the  solution  of 


~—(l+n  +  s)-0  (39) 


The  other  stability  boundary,  /  -  0,  occurs  at  n  -  3  -  s.  A  map 
of  the  linear  stability  appears  in  Figure  6  for  downflow  with 
countercurrent  ( s  <  0)  or  cocurrent  (s  >  0)  shear.  The  relative 
positions  of  the  three  boundaries  are  identical  to  those  for  the 
free-falling  film  problem.  They  intersect  at  s  -  1,  n  -  2,  which 
does  not  correspond  to  any  physical  situation  (T—  ±°°).  For 
upward  flow  (T  <  —  I,s  >  1)  the  relative  positions  of  the  three 


Figure  6.  Stability  map  for  smooth  film  in  the  presence  of 
shear  for  two-variable  case. 
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bifurcations  are  interchanged.  The  phase  plane  analysis  for 
downflow  with  countercurrent  shear  is  analogous  to  the  T  -  0 
case.  At  n  -  n„,  H  -  1  is  a  solution  for  any  w  and  the/  -  0  curve 
acquires  the  shape  of  Figure  4a.  Changing  n  modifies  the  dia¬ 
gram  in  a  similar  fashion  and  we  conclude  again  that  a  homo¬ 
clinic  orbit  must  exist  close  to  nH  with  the  same  general  shape  as 
in  Figure  5.  The  substrate  height,  H},  and  the  limiting  wave 
amplitude,  //,,  are  given  by 


Laminar  Films  with  Surface  Forces:  General  Case 

When  surface  tension  effects  arc  retained  in  the  model  a  gen¬ 
eral  form  of  Eq.  19  can  be  written  as  follows  where  w  -  H,  and 

<t>  -  tv,: 


4>i  ~ 


/(//,  w)  - 

<,«3 


g(H,  w,<(>) 
<>//3 


(43) 


-I  +  Vl  +  4 («w  -  1)(1  -  s) 
2(1  -  s) 


—  1  +  Vn!  +  12/i  —  1 
2s 


(40) 


This  analysis  shows  the  following  behavior  of  the  system  for 
downflow  with  countercurrent  shear 

1.  The  wave  characteristics  depend  on  the  shear  parameter, 
T  (or  s) 

2.  The  wave  velocity,  n,  which  approximates  the  velocity  at 

the  Hopf  bifurcation,  varies  weakly  with  T  (or  s),  taking  on  val¬ 
ues  of  1.69  at  T'-O(j-O),  1.5  at  7'-  -  —  !),and  1.75  at 

T  **  — J/«  (i  _  -3) 

3.  The  maximum  wave  amplitude  increases  markedly  with 
countercurrent  shear  since  H5  increases  as  s  becomes  more  neg¬ 
ative 

Extension  of  the  analysis  to  upflow  shows  that  the  relative 
positions  of  the  bifurcation  points.  Figures  2,  6,  are  arranged  in 
a  mirror  image  to  the  downflow  case  and  the  conclusions  are 
thus  similar  for  cocurrent  flow.  However  the  existence  of  multi¬ 
ple  values  of  the  equilibrium  film  thickness  for  any  shear  stress 
and  flow  rate  as  discussed  above  raises  these  additional  issues: 

1 .  Either  one  or  both  of  these  solutions  may  not  be  stable 

2.  The  location  of  the  turning  point  at  n  -  (3  -  T)/ 4  may  be 
outside  the  domain  where  oscillations  can  exist 


Turbulent  Film  Flow 

The  analysis  of  Eq.  19  for  laminar  film  flow  presented  above 
showed  that  C,(H  -  1 )  -  0  is  a  necessary  condition  for  instabil¬ 
ity  (see  also  Figure  6).  Note  that  Ct  is  the  coefficient  of  Hh  the 
slope  of  the  wave  in  dimensionless  form,  and  H  -  1  is  the  condi¬ 
tion  for  the  smooth  film,  h  -  h„Q  -  Qr.  For  turbulent  flow  the 
same  condition  can  be  shown  to  exist.  Thus,  the  coefficient  of  the 
dimensional  slope,  h( ,  in  Eq.  23  must  be  zero,  resulting  in  the 
following  relation  for  the  dimensionless  wave  velocity. 


where  w  -  Ht  and  <t>  -  h>,  and /(//,  w)  is  defined  in  Eq.  30.  The 
Jacobian  matrix  at  H  -  1 ,  w  -  <j>  -  0  is: 


J  - 


I 


\ 


0 

0 

QgO) 


1  o  ’ 

0  1 

/ 


(44) 


The  characteristic  equation,  \J  -  XI  |  -  0,  is 


^^(,),!  +  y±lMU^.0  (45) 


and  its  roots  are  the  eigenvalues  that  determine  the  stability. 
Applying  Routh-Hurwitz  stability  criteria  reveals  that  H  -  1  is 
stable  (X,,  X2,  Xj  <  0)  if  the  following  three  conditions  are  satis¬ 
fied 


(a)  CU  0>0 

(b)  C2(1)>0 

, ,  Q«0)  (I  +  T)2c,(\)CA\h2  „ 

(c)  - - - >  0 

«i  A 


(46) 


Condition  (a)  is  the  exchange  of  stability,  described  earlier,  and 
it  implies  n  <  3  -  j.  Condition  (b)  is  also  unchanged,  requiring 
n  >  (3  +  s)j 2.  Hopf  bifurcation  occurs  when  condition  (c)  is 
violated  and  the  bifurcation  points  depends  now  on  the  ratio  of 
the  viscous  to  surface  tension  terms.  It  reduces  either  to  the  pre¬ 
vious  Hopf  condition  (C,  -  0)  when  <2/<,  —  °°  or  to  condition  (a) 
when  e2/r,  — *  0.  Substituting  the  relations  for  the  Q(H),  Eq.  20, 
we  find  that  Hopf  bifurcation  (n  -  nH)  occurs  at 

0  (1  +  T)*,}Rei,} 


-it2  +  2/(a)n  -/(a)  -  0  /(«)  -  (41) 

a(2  +  a) 

or 

«-/(«)  +  V/(/-  0  (42) 

This  result  is  independent  of  the  particular  form  of  the  correla¬ 
tion  used  for  the  wall  stress.  At  high  flow  rates,  as  the  flow 
becomes  more  turbulent  a  increases  and  n  approaches  1 .0.  Thus 
for  a  velocity  distribution  that  follows  a  ‘/th  power  law,  n  - 
1 .14.  It  should  be  noted  that  such  low  values  of  the  wave  velocity 
cannot  be  predicted  from  any  laminar  distribution  but  have  been 
observed  experimentally  for  high  flow  rates,  Eq.  12. 


- (3~J  -  n)'° _  (47) 

3  s\ 

n  -  --  -\[-2n2  +  12 (n  -  l)2  -  s(l  +  n  +  j)] 

The  nature  of  the  steady  state  is  shown  in  the  (/3,  n)  plane  of 
Figure  7  by  denoting  the  sign  of  the  three  eigenvalues.  Limit 
cycles  may  exist  in  the  range  0  +  s)/2  <  n  <  nH  and  the  order  of 
the  various  points  is  unaffected  by  changing  T  for  downward 
flow. 

Analysis  of  the  trajectories  in  the  three-dimensionai  space 
(//,  w,  4>)  is  more  intricate  than  the  two-dimensional  version. 
The  structure  of  the  analysis  follows.  We  investigate  the  shape 
of  the  g  -  0  surface,  Eq.  43,  in  the  phase  space  to  show  that  it 
attracts  the  motion  and  ask:  is  the  motion  along  it  stable?  It 
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Figure  7.  Stability  map  for  general  case. 


turns  out  that  it  is  unstable  in  certain  sections.  We  find  the  sta¬ 
ble  and  unstable  manifold  and  construct  the  trajectory  in  space. 
The  steady  state  H  -  I  may  be  the  origin  of  the  homoclinic 
orbit.  It  has  a  stable  one-dimensional  (i.e.,  a  line)  manifold  and 
an  unstable  two-dimensional  (surface)  manifold  along  which 
the  motion  spirals.  These  spirals  yield  the  capillary  wave  and  are 
evident  only  with  relatively  small  Re  (t,/e 2  large).  After  it  spi¬ 
rals  out,  the  trajectory  travels  in  the  vicinity  of  Hy  so  that  the 
dimensionless  value  of  the  film  thickness  is  still  a  good  approxi¬ 
mation  for  the  substrate  thickness. 

The  slopes  of  the  trajectories  can  be  found  readily  from  Eq.  43 
to  be: 

**  _  d±  _ 

dH  tlH’w  dw  t,//30  v  1 

Note  that  t,  — *  0  except  for  the  lowest  Reynolds  numbers. 
Therefore  the  slopes  of  the  trajectories  in  both  the  0  -  H  and 
<(>  -  w  planes  are  steep  everywhere  in  the  phase  space  except 
where  g— »  0  or  when  the  first  and  second  derivatives  of  the  wave 
profile,  k>  and  0 ,  are  large.  Data  show  that  these  derivatives  are 
everywhere  very  small  indeed.  Therefore  we  now  explore  the 
region  around  g  -  0. 

From  the  definition  for  g  in  Eq.  43,  when  g  -  0 


f{H,  w) 
<jCj(W) 


(49) 


1  H 


Figure  8.  Family  of  f-0  curves  in  w-H  plane. 

can  be  described  by  the  relations 

H,-  w  h-,-0-/A3C3(//)  (50) 

This  is  identically  the  two-dimensional  model  developed  ear¬ 
lier,  but  the  question  of  stability  at  any  position  on  the  surface 
must  still  be  explored. 

Suppose  He(l),  w„(l),  0„(/)  is  the  solution  to  Eq.  43  along  g  - 
0  with  perturbations  from  this  motion  H  -  //(/)  -  //„(/),  w  - 
h '(/)  -  w,(l)  and  0  -  0(/)  -  0O(/)  described  by 

H,-w  tv,  -  0  -  g„H  +  gj.  +  g,0  (51) 


Assume  that  perturbations  in  H  are  small  and  grow  slowly  sr 
that  g„  -  dg/dH  can  be  averaged,  and  (gH),  {gj,  and  (gt) 
replace  the  original  time-dependent  gH,  gw,  gf.  The  motion  is  sta¬ 
ble  if  the  trivial  solution  of  Eq.  51,  H  -  0  -  iv  -  0  is  stable,  that 
is,  if  the  eigenvalues  are  al!  negative.  We  do  a  stability  analysis 
similar  to  that  leading  to  Eq.  46  to  find  the  stability  condition 


<g*> 

r,//1 


<g*><gw>  „ 

-  c  0 


(52) 


For -  0  this  implies  that /-  0,  a  situation  that  was  explored 
earlier.  This  can  be  pictured  as  a  cylindrical  surface  in  the  phase 
space  perpendicular  to  the  w  -  H  plane.  Figure  8  shows  f-0 
curves  for  a  sequence  of  values  of  n  as  determined  from  Eq.  37. 
The  case  for  rt  -  1 .69  appeared  earlier  in  Figure  4a  where  it  was 
shown  that  the  crossing  of  the  abscissa  m:  rks  the  location  of  H}, 
the  substrate  height.  Note  that  when  n  -  3  the  value  of  H, 
becomes  identical  with  the  smooth  film  thickness,  a  condition  in 
agreement  with  the  experiment. 

Figure  9  pictures  the  cylindrical  surface,  g  -  0,  in  the  phase 
space  0,  w,  H  for  «3  -  0.  When  t small  but  not  zero  the  surface 
is  tilted  somewhat.  All  trajectories  starting  frem  any  point  in  the 
phase  space  not  on  g  -  0  rapidly  approach  this  plane  since 
drp/SH  and  d<t>/dw  are  large.  This  is  evident  from  Eq.  48  for  the 
condition  r  — *  0.  Once  the  trajectory  reaches  the  g  -  0  surface  it 


Figure  9.  g  -  O  surface  In  phase  plane. 
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An  analysis  of  this  equation  shows  that  the  transition  occurs  at 
H  -  1-  The  section  H  <  1  is  stable  and  is  the  stable  manifold 
near  H  -  1 .  The  motion  around  H  >  1  is  unstable  and  the  trajec¬ 
tory  oscillates  around  the  line /-  0. 

We  now  describe  the  motion  in  the  vicinity  of  H  -  1  with  n  ~ 
nH.  The  stable  manifold  is  the  portion  of  the  curve /  -  0,  denoted 
act  in  Figure  10,  that  connects  the  Ht  and  H  -  1  solutions, 
denoted  ABC.  Its  trajectory  flows  toward  the  steady  state.  At 
n  -  nH  (Hopf  point)  there  is  a  pa:r  of  imaginary  eigenvalues  and 
one  negative  eigenvalue: 


X|  — - C3(l),  Xjj  —  ±  i 

<i 


(54) 


The  oscillations  develop  or  decay  with  a  period  of 
r,(l))1/2.  Since  C,(l)  ~  short  periods  of  0(  \f7x)  are 

expected  when  the  surface  tension  is  important,  while  when  sur¬ 
face  tension  effects  are  negligible  the  period  is  of  0(  V^).  The 
eigenvectors  that  correspond  to  these  imaginary  eigenvalues 
were  calculated  to  be  <t>  ■*  -  1 )  and  <t>  -  ±r'6w;  they  define 

the  unstable  manifold.  We  note  however  that  62  — ►  oo  as  <,  — ►  Oor 
— ►  0.  Thus  52  (W  -  1 )  is  finite  only  for  H  -  1 .  The  plane  H  -  1 
is  thus  c  good  approximation  for  the  unstable  manifold.  For  n 
somewhat  smaller  than  n,„  X,  does  not  change  significantly, 
while  Xjj  -  ±  i&.  The  frequency  is  of  order  1/2  or  «j 1/2  as 

before  and  'F  is  the  amplitude  growth  rate,  which  is  small  for 
high  surface  tension  and  large  at  the  other  extreme.  For 
extremely  small  values  <  f  (rtH  -  n)  a  limit  cycle  exists  and  it  lies 
close  to  the  H  -  1  plane.  Further  increase  in  (nH  -  n)  will  break 
the  cycle  and  the  trajectory  is  detached  from  the  plane  and 
moves  fast  toward  the  stable  manifold  (/  -  0).  The  resulting 
motion  in  Figure  10  shows  a  closed  curve  composed  of  a  spiral¬ 
ing  motion  on  H  -  1  followed  by  a  loop  that  travels  around  Hy 
In  the  time  domain  we  find  small  capillary  waves  whose  number 


Table  I.  Wave  Characteristics  of  Free-Falling  Film  (T  -  0) 


n  Substrate 

Qr/t’o  c - - 


Re 

0 

r.irn 

cm/s 

cm/s 

Meas. 

Theory 

Meas. 

Theory 

78 

1.18 

0.28 

28.2 

68 

2.58 

2.46 

0.8 

0.75 

192 

0.26 

0.39 

50.3 

105 

2.14 

2.10 

0.72 

0.72 

387 

0.08 

0.50 

76.8 

120 

1.49 

1.77 

0.66 

0.62 

778 

0.025 

0.63 

127.0 

140 

!.I0 

1.70 

— 

47 

Tabic  2.  Measured  Conditions  Along  Flooding  Curve 


t,  c  h , 


Re 

N/m2 

T 

-  5 

cm/s 

n 

cm 

15 

3  29 

0.9778 

44  1 

35.3 

11.4 

0017 

25 

3.28 

0.9662 

28.5 

43.8 

905 

0  0215 

43 

3.32 

0.9237 

12  1 

500 

6.29 

0024 

63 

3.49 

0.9041 

9  4 

58  2 

5  31 

0026 

78 

3.25 

0.8906 

8,1 

61.6 

3  60 

C.029 

109 

3.59 

0.8722 

6.8 

56  0 

3.20 

0  030 

189 

3.80 

0.8309 

4.9 

86  9 

3.14 

0037 

298 

3.98 

0.7702 

3,35 

96.36 

2  56 

0.C42 

varies  like  <,/cj,  Figure  5.  The  trajectory  disappears  when  it  hits 
H3,  which  is  still  a  good  approximation  for  substrate  thickness. 

This  analysis  of  the  general  third-order  equation  showed 
again  that  n„  and  H}  arc  good  approximations  for  wave  velocity 
and  substrate  thickne^.  respectively,  but  that  the  wave  shape 
differs  from  that  of  the  second-order  model. 

Comparison  with  Experiments 

Experimental  data  on  wave  structure  on  falling  liqu'd  films 
with  counterflow  of  gas  were  recently  presented  by  Zabaras  and 
Dukler  (1988).  These  data  are  used  here  to  evaluate  this 
model. 

Free-falling  film 

Equation  47  with  s  -  0  can  be  used  to  determine  nH,  the 
dimensionless  wave  velocity  at  Hopf  bifurcation,  and  this 
can  be  compared  with  the  experimental  values  to  test  the  conclu¬ 
sion  of  the  phase  space  study  that  this  must  be  close  to  the  veloc¬ 
ity  of  the  homoclinic  orbit  and  thus  to  experimental  data.  Simi¬ 
larly,  Eq.  40  provides  a  value  of  H3,  the  minimum  possible  film 
thickness  computed  from  the  theory.  This  can  be  useo  to  com¬ 
pute  the  substrate  thickness,  which  can  be  compared  with  the 
data.  These  comparisons  with  the  data  of  Zabaras  and  Dukler 
(1988)  appear  in  Table  1;  surprisingly  good  agreement  is  dis¬ 
played  considering  the  difficulty  in  experimental' v  determining 
the  substrate  thickness.  It  should  be  noted  that  the  Reynolds 


Figure  11.  Theory-experiment  comparison  tor  wave  ve¬ 
locity  along  flooding  curve. 
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Figure  12.  Theory-experiment  comparison  for  substrate 
thickness  along  flooding  curve. 

o  Theory;  •  Experiment 


/.  -  length  scale  h,Rrj(\  *  1 ) 

,V  -  wave  number 

n  -  dimensionless  wave  lelncity  chJQ, 
n„ .  n,  -  n  of  the  Hopf  and  homoclimc  saddle  loop  bifurcations 
P  -  pressure 

Q.  Qf  -  local  and  feed  liquid  volumetric  flow  rate 
Re  -  Reynolds  lumber.  Qjv 
s  -  modified  shear  parameter  77(  I  ^  T) 

T  -  dimensionless  interfacial  shear,  3  r,/2  pg  h„ 

/  -  time 

Uc  -  linear  gas  velocity 

u,  *>  -  liquid  velocity  parallel,  perpendicular  to  wall 
w  -  dff/dl 

W  -  surface  tension  parameter.  Eq.  18 
ff'e.  B7  -  gas,  liquid  mass  flow  rate 

x,  y  -  coordinates  in  direction  parallel,  perpendicular  to  wall 


numbers  in  this  work  are  defined  as  one-fourih  those  in  the 
Zabaras  and  Dukler  paper. 

At  Reynolds  number  above  about  300  turbulent  flow  is 
thought  to  set  in,  the  shape  of  the  velocity  profile  changes,  and 
the  wave  velocity  would  be  expected  to  decrease,  as  discussed 
earlier. 

Film  with  countershear 

A  severe  test  of  the  model  is  a  comparison  with  data  taken 
under  the  condition  of  interfacial  shear  due  to  counterflow  of 
as.  Table  2  lists  some  original  source  data  taken  from  the  thesis 
of  Zabaras  (1985)  during  experiments  along  the  flooding  curve. 
Note  that  in  accordance  with  theory  the  value  of  the  dimension¬ 
less  interfacial  shear,  T,  is  never  less  than  —  1 .  However  at  low 
liquid  rates  it  approaches  it  closely.  Figure  11  compares  the 
experimental  values  of  n  with  nH  calculated  from  the  simple 
form  of  the  equation  in  the  limit  /3  — -  0.  Eq.  38,  and  from  the 
general  analysis,  Eq.  47.  The  comparison  is  most  encouraging, 
especially  considering  the  fact  that  when  T  — *  -1,  the  com¬ 
puted  value  of  s  -  T/(  1  +  T)  used  in  both  of  these  equations  is 
extremely  sensitive  to  small  errors  in  T. 

Not  quite  so  satisfactory  a  comparison  is  shown  in  Figure  1 2. 
Here  the  dimensionless  thickness  calculated  from  Eq.  40  is  com¬ 
pared  with  the  dimensionless  substrate  thickness.  The  experi¬ 
mental  values  of  j  and  n  were  used  as  input  to  these  calculations. 
The  error  bars  on  the  result  due  to  the  sensitivity  of  s  to  errors  in 
T  at  the  low  flow  rate  is  indicated.  However  it  is  still  clear  that 
the  theory  is  not  quite  accurate  enough  in  this  region.  We  sus¬ 
pect  this  deviation  is  due  to  the  assumption  that  the  interfacial 
shear  stress  is  independent  of  position  along  the  interface  and 
thus  is  the  same  over  the  substrate  as  it  is  over  the  wave  itself. 
Experiments  suggest  that  this  is  not  the  case  and  we  are  now 
exploring  revised  models  which  incorporate  shear  stress  varia¬ 
tion  along  the  wave. 
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Notation 

C,  -  functions,  Eq.  20 
c  -  dimensional  wave  velocity 
D  -  pipe  diameter 

/*/,  -  parameter  in  r,  correlations,  Eq.  58 
/,  g  -  functions,  Eqs.  30, 43 
h.  H  -  dimensional,  dimensionless  film  thickness 
h „  h„  -  smooth  film  solution  with,  without  shear 
J  -  Jacobian  matrix 
/  -  dimcnft.'inless  coordinate  f/L 


Greek  letters 

a  -  exponent  in  power  law  profile.  Eq.  21 

d  -  ratio  of  time  scales,  Eq  47 

e  -  kinematic  viscosity 

p  -  density 

t,  -  shear  at  interface 

a  -  surface  tension 

<1  -  dimensionless  coordinate,  y/h 

it  -  viscosity 

0.  (  -  lime,  distance  in  a  moving  frame 
<]  -  dimensionless  time  scales  in  terms  Huh  H„.  defined  fol¬ 
lowing  Eq.  43 
<j>  -  d'H/dP 

X  -  eigenvalues  of  Jacobian  mail  x 

Subscripts 

x,  y.  H,  etc.  -  derivative  with  respect  to  denoted  quantities 
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Abstract — Mass  transfer  from  a  solid  boundary  into  a  thin  wavy  film  was  studied  for  a  wide  range  of 
laminar  flow  rates  with  a  novel  experimental  technique.  Speculations  that  large  waves  control  the  transport 
process  through  a  convection  mechanism  are  validated  through  an  examination  of  the  time  variation  of  the 
film  thickness  and  solute  concentration.  Statistical  analyses  of  the  data  demonstrate  that  the  occurrence  of 
large  waves  with  higher  local  flow  rates  coincides  with  an  increased  mass  transfer  both  locally  and  globally. 
Accordingly,  the  relative  importance  of  small  waves  or  ripples  is  shown  to  decrease  rapidly  with  increasing 
flow  rate. 


INTRODUCTION 

As  a  liquid  film  flows  down  a  solid  surface,  a  wavy 
interface  develops.  For  flow  rates  of  industrial  inter¬ 
est,  the  free  surface  quickly  becomes  covered  with  a 
complex  array  of  waves  whose  amplitudes  vary 
greatly  about  the  mean  film  thickness.  This  gravity- 
driven  behavior  is  observed  even  in  the  absence  of 
interfacial  shear  stresses  induced  by  adjacent  gas  flow. 
Figure  1  shows  a  sample  time  tracing  of  the  interface 
of  a  fully  developed  laminar  film  falling  freely  down  a 
vertical  tube.  Large  waves,  having  amplitudes  from  2 
to  5  times  the  mean  thickness,  move  over  the  thin 
substrate  at  velocities  up  to  several  times  the  mean 
and  carry  a  large  fraction  of  the  total  liquid  flow.  The 
slope  of  these  large  waves  seldom  exceeds  5%,  yet 
numerical  studies  of  the  flow  fields  (Wasden  and 
Dukler,  1989)  show  that  significant  normal  velocities 
exist  near  the  gas-liquid  and  wall  interfaces. 

Mass  transfer  in  liquids  is  a  slow  process  compared 
to  heat  transfer,  suggesting  that  mass  transfer  rates 
respond  more  dramatically  than  heat  transfer  rates  to 
wave-induced  fluctuations  in  the  otherwise  parallel 
flow  (cf.  Seban  and  Faghri,  1978;  Henstock  and 
Hanratty,  1979).  Many  analytical  studies  of  transport 
through  flat  or  slightly  rippled  films  exist,  all  treating 
the  interface  as  a  regular  or  periodic  surface.  For 
selected  large  waves,  the  convective  fluxes  due  to 
normal  velocities  near  the  interfaces  dominate  the 
transport  and  partially  explain  why  the  previous  ana¬ 
lyses  underptedicted  the  overall  mass  transfer  rates 
(Wasden  and  Dukler,  1990). 

Previous  experimental  studies  of  mass  transfer  into 
falling  films  have  demonstrated  that  transfer  rates  in 
wavy  films  are  greater  than  when  waves  are  sup¬ 
pressed,  without  providing  an  insight  into  the  mech¬ 
anisms  of  enhancement.  Gas  absorption  experiments 
are  typically  conducted  by  contacting  the  entire  film 
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surface  with  an  excess  of  solute  and  measuring  the 
bulk  outlet  concentration  (Emmert  and  Pigford, 
1954).  This  practice  yields  data  averaged  over  all 
waves  on  the  surface  and  makes  it  impossible  to 
explore  the  effects  of  wave  characteristics.  An  altern¬ 
ative  arrangement  to  study  wave  effects  on  mass 
transfer  is  illustrated  in  Fig.  2.  Solvent  flows  as  a  wavy 
film  over  a  contaminated  wall  section,  which  may  be 
located  any  distance  below  the  solvent  entrance. 
Large  waves  evolve  slowly  on  the  fully  developed 
films;  so  if  the  film  thickness  and  the  average  contam¬ 
inant  concentration  are  measured  directly  after  the 
film  contacts  this  section,  the  effects  of  wave  shape, 
size,  velocity,  etc.,  on  the  mass  transferred  into  the  film 
should  be  plain.  According  to  the  criterion  of  Arts 
(1956),  Taylor  dispersion  in  this  system  is  negligible. 
This  configuration  is  similar  to  the  one  studied  by 
Stirba  and  Hurt  (1955)  and  offers  several  experimental 
advantages  for  studying  the  effect  of  waves.  First, 
while  the  presence  of  waves  causes  interfacial  transfer 
rates  for  gas  absorption  to  be  2-4  times  greater  than 
those  observed  for  smooth  films  (Henstock  and 
Hanratty,  1979),  mass  transfer  from  the  wall  appears 
to  occur  5-10  times  faster  than  in  smooth  films  (Stirba 
and  Hurt,  1955).  Second,  it  is  possible  to  maintain 
precise  control  over  the  amount  of  contaminant 
entering  the  solvent,  whereas  gas  absorption  opera¬ 
tions  are  limited  to  providing  contaminant  in  excess. 
Third,  the  effects  of  wave  amplitude  and  velocity  may 
easily  be  studied  by  varying  the  distance  between  the 
film  entry  and  the  contaminated  section. 

Previous  experimental  studies  of  this  type,  which 
used  the  dissolution  of  a  solid  solute  surface  by  the 
film  (Stirba  and  Hut,  1955;  Kramers  and  Kreyger, 
1956;  Oliver  and  Atherinos,  1968),  measured  mass 
transfer  rates  either  by  weighing  the  flow  surfaces 
before  and  after  an  experiment  or  by  measuring  the 
bulk  concentration  of  solute  in  the  fluid  leaving  the 
apparatus.  Such  procedures  made  it  impractical  to 
explore  the  effects  of  wave  characteristics.  Simultan¬ 
eous  measurements  of  mass  concentration  and  film 
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Fig.  1.  Film  thickness  time  trace,  Re  =  880. 


Fig.  2.  Mass  transfer  from  the  solid  boundary  into  a  passing 
wave. 


structure  have  not  been  reported  for  transfer  from  the 
wall  into  a  film,  but  Stainthorp  and  Wild  (1967) 
measured  the  film  thickness  and  the  concentration 
averaged  over  the  film  thickness  in  a  thin-film 
liquid-liquid  extraction  column.  Measurements  were 
taken  using  two  photocells  located  38  mm  apart  in  the 
flow  direction  and  an  absorbance  technique  was  used 
to  relate  photocell  readings  to  film  thickness  and 
concentration. 


The  analysis  of  mass  transfer  from  the  wall  io  a 
liquid  is  represented  by  a  substantial  body  of  literat¬ 
ure  due  to  the  interest  in  using  electrochemical  probes 
with  a  redox  reaction  of  measure  wall  shear  stress  in 
high  Schmidt  number  fluids  (see  the  review  by 
Hanratty  and  Campbell,  1983).  Under  these  condi¬ 
tions,  the  concentration  boundary  layer  near  the  wall 
is  exceedingly  thin  and  the  velocity  profile  is  con¬ 
sidered  linear.  For  quasi-steady  conditions  a  closed- 
form  relation  between  mass  transfer  rale  and  wall 
shear  stress  can  be  derived.  However,  this  method  is  of 
little  use  when  the  resistance  to  mass  transfer  extends 
through  the  film  as  in  large  waves. 

Experimental  studies  have  examined  the  average 
effects  of  waves  on  transport  and  numerical  ap¬ 
proaches  have  explored  transport  processes  in  isol¬ 
ated  waves.  Both  avenues  show  the  importance  of 
wave  characteristics  and  suggest  the  need  for  detailed 
simultaneous  local  measurements  of  film  structure 
and  concentration.  The  work  presented  here  includes 
the  description  of  a  novel  method  to  obtain  these 
measurements  and  an  analysis  of  the  laminar  film  flow 
results. 


experimental  apparatus,  technique  and 

PROCEDURE 
Flow  loop  and  measuring  station 

The  flow  loop  is  illustrated  in  Fig.  3.  Solvent  enters 
the  top  of  the  vertical  column  through  a  stainless  steel 
porous  sinter  which  distributes  the  inlet  flow 
smoothly  and  uniformly  around  the  tube  periphery. 
Solute  is  pumped  through  the  pores  of  a  very  smooth 
plastic  cylindrical  insert  located  directly  above  the 
measuring  station.  The  column  is  constructed  of  a 
number  of  precisely  machined  segments  of  Plexiglas 
tubing  with  50.8  +  0.2  mm  ID.  The  flow  development 
distance  is  adjusted  by  varying  the  number  of  seg¬ 
ments  between  the  solvent  feed  and  solute  input  sec¬ 
tions.  Solvent  (an  aqueous  salt  solution)  and  solute 
(an  aqueous  salt  and  dye  solution)  were  pumped  into 
the  system  using  magnetically  coupled  gear  pumps. 
The  pump  drive  units  contain  an  internal  feedback 
circuit  to  maintain  constant  flow  rates  with  minimal 
fluid  heating.  Details  of  the  solvent  feed  section, 
pumps  and  flow  loop  construction  arc  found  else¬ 
where  (Wasden,  1989). 

The  solute  feed  section  was  designed  to  uniformly 
distribute  a  very  low  flow  rate  of  contaminant  both 
axially  and  circumferentially  so  that  the  section  ap¬ 
proximates  a  contaminated  solid  wall.  This  arrange¬ 
ment  was  developed  to  avoid  previously  reported 
difficulties  in  wetting  a  solid  contaminant  surface 
(Stirba  and  Hurt,  1955;  Kramers  and  Kreyger,  1956; 
Oliver  and  Atherinos,  1968).  As  shown  in  Fig.  4,  the 
section  is  made  of  an  extruded  piece  of  porous  poly¬ 
ethylene  whose  surfaces  are  very  smooth.  Surface  pore 
sizes  average  10  pm  and  solute  is  constantly  pumped 
at  very  low  rates  from  the  annulus  through  the  pores 
into  the  solvent  stream.  The  porous  plastic  ID  was 
slightly  larger  than  the  Plexiglas  tube  ID,  so  the 
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Fig.  3.  Experimental  flow  loop. 


Fig.  4.  Contaminant  feed  section  of  flow  loop. 


flanges  connecting  the  280  mm  long  porous  section  to 
the  column  were  smoothly  tapered  at  a  2.9°  angle. 
Uniform  wetting  of  the  surface  was  achieved  by 
blocking  a  large  portion  of  the  back  side  of  the  porous 
plastic.  Fluid  enters  the  solute  section  through  the 
uncovered  portion  of  the  outer  surface.  Once  past  this 
barrier,  hydrostatic  and  pump  pressure  equalize  the 
fluid  distribution  in  the  highly  porous  central  region 
and  force  it  to  the  inner  wall,  where  it  uniformly 
saturates  the  small  pores.  By  adjusting  the  area  closed 
to  flow  on  the  outside  of  the  porous  section,  fluid  was 


observed  to  be  distributed  evenly  both  axially  and 
circumferentially  over  the  lower  255  mm  of  the  device. 

The  measuring  station  (Fig.  5)  included  provisions 
for  measuring  the  film  thickness  at  two  axial  locations 
and  the  average  concentration  at  the  upper  film  thick¬ 
ness  location.  The  device  was  constructed  from  a  piece 
of  Plexiglas  bored  to  50.8  mm  in  which  two  sets  of 
parallel  76  nm  diameter  Ni-Pt  wires  axially  separated 
by  25.0  mm  were  installed  to  measure  the  film  thick¬ 
ness.  A  set  of  coilinear  quartz  rods  with  highly 
polished  ends  was  used  to  transmit  a  laser  light  beam 
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Fig.  5.  Wave  structure  and  contaminant-measuring  section. 


through  the  contaminated  film.  The  center  of  the 
receiving  rod  was  mounted  flush  with  the  wall,  be¬ 
tween  and  1.7  mm  below  the  top  wires.  The  transmit¬ 
ting  rod  penetrated  the  film  on  the  opposite  side  of  the 
tube.  One  end  of  this  rod  was  positioned  so  that  it  was 
12.5  mm  from  the  receiving  rod  end.  In  order  to 
minimize  reflective  losses  from  the  liquid  surface, 
12.5  mm  diameter  parabolic  mirror  with  a  12.5  mm 
focal  length  was  placed  over  the  rod  end.  A  pair  of 
100  pm  holes  drilled  above  the  central  mounting  hole 
permitted  the  parallel  conductance  wires  to  pass 
through.  A  set  of  “O”  rings  was  placed  on  the  outside 
of  the  penetrating  rod  to  keep  liquid  from  migrating 
from  the  film  onto  the  rod  end  and  obscuring  the  light 
path. 

Measurement  techniques  and  instrumentation 
Two  experimental  techniques  were  combined  in 
this  study.  A  high-frequency  circuit  was  used  to  relate 
the  conductance  of  the  liquid  between  the  parallel 
wires  to  the  film  thickness.  The  resolution  of  the 
method  is  of  the  order  of  15  pm  and  the  circuit  had  a 
time  lag  of  less  than  100  ps  (Brown  et  ai,  1978).  The 
average  concentration  was  determined  by  measuring 
the  attenuation  of  a  laser  light  beam  passing  through 
the  contaminated  film  while  simultaneously  meas¬ 
uring  the  film  thickness  with  the  parallel-wire  con¬ 
ductance  method.  The  relation  between  optical  at¬ 


tenuation  and  solute  concentration  invokes  Beer’s  law 
for  monochromatic  radiation  (Knowles  and  Burgess, 
1984).  For  a  single  contaminant, 

A  =  eh<C>  (1) 

where  A  is  the  optical  absorption  of  the  contaminated 
solution,  c,  the  molar  absorptivity  of  the  contaminant 
at  the  particular  wavelength,  h,  the  path  length  of  the 
radiation  and  <C>,  the  average  concentration  of  the 
contaminant  over  the  path  length.  The  absorption  can 
be  related  to  the  incident  and  transmitted  radiation 
through 

A  =  In  (Jo/I)  «  In  (1/TJ  (2) 

where  J0  is  the  intensity  of  the  incident  beam,  /,  the 
attenuated  intensity  and  Tt,  the  transmittance  of  the 
medium  for  the  wavelength  of  interest.  An  optically 
dense  dye  solution  was  used  as  the  contaminant.  By 
using  a  transparent  aqueous  salt  solution  as  the 
solvent,  the  attenuation  of  light  indicates  only  the 
amount  of  solute  in  the  optical  path. 

A  schematic  of  the  optical  measurement  system  is 
shown  in  Fig.  6.  The  light  source  was  a  5  mW  He-Ne 
laser  (632.8  nm)  concentrated  into  a  0.8  mm  diameter 
beam  and  rated  at  less  than  0.1%  RMS  noise.  A  6" 
diameter  integrating  sphere  was  used  to  capture  light 
leaving  the  quartz  rod  and  provide  a  uniform  lumi¬ 
nous  flux  on  the  face  of  a  photodetector  mounted  90° 
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Fig.  6.  Optical  measurement  system. 


F 


Mass  transfer  from  a  wall  into  a  wavy  failing  film 


4327 


I 

I 

I 

I 


I’ 

sriitlted 

uninant 


off-axis  of  the  incoming  beam.  This  arrangement  in¬ 
sured  that  the  output  of  the  photodetector  was  not 
influenced  by  the  alignment  of  the  photodetector  face 
with  respect  to  the  quartz  rod.  The  photodetector 
current  indicates  the  light  transmitted  and  is  ex¬ 
pressed  as 

Tx  =  exp  ( —  A)  =  exp  ( —  eh(C))  m 

l-eh(C)  ■■■  (3) 

for  small  values  of  absorption.  The  value  of  <C)  was 
determined  once  h  was  obtained  from  the  parallel  wire 
probes.  Dye  concentrations  were  chosen  so  that  the 
transmittance  was  always  greater  than  95%  and  the 
linear  approximation  [eq.  (3)]  was  accurate  to  within 
0.1%.  Analog  data  collected  from  the  film  thickness 
probes  and  photodetector  circuitry  were  digitized  and 
subsequently  processed.  All  signals  were  digitized 
at  a  frequency  of  1  kHz  per  channel  with  a  12-bit 
analog-digital  converter  installed  in  a  micro¬ 
computer. 

Solvent  and  solute  properties 
The  solvent  was  a  0.34  M  aqueous  solution  of  NaCl 
with  0.04  gl~'  dye  The  dye  used  was  toluidine  blue 
(Basic  Blue  17)  which  has  a  peak  in  its  absorption 
spectra  at  628  nm.  The  molar  absorptivity  of  the  dye 
(c)  at  630 nm  is  quoted  as  5000 M "’em'1.  This 
highly  absorbing  dye  was  chosen  so  that  a  very  small 
concentration  could  be  detected.  An  aqueous  solution 
of  5  g  1  * 1  dye  and  0.34  M  NaCl  was  used  as  the  solute. 
The  solute  dye  concentration  was  kept  low  to  insure 
that  the  physical  properties  of  the  solution  were  not 
appreciably  altered  by  its  presence.  At  22°C  the 
solvent  had  the  following  measured  properties: 
density,  1010  kgm~3,  absolute  viscosity,  9.845  x  10~4 
kgm'ls'‘  and  surface  tension  coefficient,  74.2 
x  10~3  Nm'1.  The  diffusion  coefficient  for  the  con¬ 
centrated  organic  dye  diffusing  into  a  slightly  contam¬ 


inated  solvent  was  estimated  to  be  2.8  x  10  10  m3s~  ' 
(Cussler,  1984).  Over  the  course  of  the  experiments  no 
deviations  greater  than  0.2'C  were  noted  and  fluid 
properties  were  considered  constant. 

Calibration ,  verification  of  measurement  techniques  and 
procedure 

Errors  in  film  thickness  measurements  were  estim¬ 
ated  to  be  no  more  than  1%.  The  error  in  the  concen¬ 
tration  averaged  over  the  film  thickness  was  shown  to 
be  less  than  2.5%.  Environmental  noise  in  the  optical 
signal  was  eliminated  by  covering  the  entire  apparatus 
in  thick  black  plastic  sheeting  and  losses  in  the  optical 
arrangement  were  accounted  for  in  the  calibration 
procedure.  The  accuracy  of  the  method  can  be  con¬ 
firmed  by  examining  data  obtained  for  the  flow  of 
pure  solvent  with  a  uniform  dye  concentration.  The 
optical  signal  is  converted  to  the  film  thickness  for  a 
known  dye  concentration  using  eq.  (3).  Figure  7 
shows  a  typical  film  thickness  trace  obtained  using  the 
conductance  method  at  a  Reynolds  number  of  750. 
Optical  measurements  were  made  simultaneously  and 
the  relative  deviation  between  conductance  and  op¬ 
tical  measurements  are  shown  in  the  figure.  The  agree¬ 
ment  proves  that  the  sensitivity  of  the  optical  tech¬ 
nique  is  sufficient  to  resolve  deviations  in  local  dye 
concentrations  of  the  order  of  1  %  of  the  solvent  or 
0.008%  of  the  solute  concentration.  The  agreement 
also  shows  that  the  frequency  responses  of  the  instru¬ 
ments  are  comparable  and,  thus,  can  measure  small 
waves. 

Each  experimental  run  was  preceded  by  calibration 
of  the  system,  the  procedures  for  which  are  described 
elsewhere  (Wasden,  1989).  Solute  flow  rates  of  roughly 
2%  of  the  solvent  flow  rate  were  found  to  provide 
adequate  contaminant  without  disrupting  the  flow. 
Data  were  collected  for  solvent  flow  rates  correspond¬ 
ing  to  Reynolds  numbers  between  100  and  1000.  Each 
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Fig.  7.  Film  thickness  and  deviations  between  optical-  and  conductance-based  measurements.  Re  ~  750. 
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solvent-solute  flow  rate  condition  was  sampled  for 
90s,  allowing  500-1000  large  waves  to  pass  the 
measuring  station. 

data  analysis  and  discussion  of  results 
Time  series  of  film  thickness  and  the  average  dye 
concentration  were  examined  to  test  the  speculation 
that  large  waves  dominate  the  transport  dynamics. 
Since  dye  concentration  was  measured  1.7  mm  below 
the  film  thickness,  cross  correlations  of  the  film  thick¬ 
ness  signals  were  computed  and  used  to  extract  a  time 
delay  for  time-shifting  the  concentration  signal  to 
approximate  simultaneous  measurements.  The  degree 


of  mixing  was  inferred  from  a  solute  distribution 
coefficient 

I  ii 

<P(l)  =  -----  I  C(y,t)dy.  (4) 

c«,n  "n  J,  =  o 

The  coefficient  would  be  identical  in  shape  to  the  film 
thickness  for  a  perfectly  mixed  film.  Conversely,  if  no 
mixing  occurs  (as  in  a  flat  film),  4>(t)  is  constant. 

Typical  time-dependent  measurements  of  the  film 
thickness  and  the  solute  distribution  coefficient,  <£, 
appear  in  Figs  8  and  9  for  Reynolds  numbers  of  230, 
where  the  large  waves  are  only  weakly  developed,  and 
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Fig.  8.  Film  thickness  and  solute  distribution  coefficient  for  a  film  development  length  of  2.14  m,  Re  =  230. 
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Fig.  9,  Film  thickness  and  solute  distribution  coefficient  for  a  film  development  length  of  2.14  m.  Re  =  950. 
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950,  where  they  dominate  the  interface.  Film  thick¬ 
nesses  are  normalized  by  hx.  the  film  thickness  com¬ 
puted  for  a  laminar  film  without  waves.  At  low 
Reynolds  numbers  the  ratio  of  the  wave  peak  to 
substrate  thickness  is  as  much  as  four  but  at  Re  =  950 
wave  peak  to  substrate  thickness  ratios  of  the  large 
waves  range  from  four  to  eight.  At  the  higher  flow  rate 
the  waves  are  asymmetrical  in  shape.  Wasden  and 
Dukler  (1989)  computed  the  flow  field  for  such  waves 
and  showed  that  for  large  asymmetric  waves  there 
exist  significant  velocities  normal  to  the  wall  at  the 
front  and  back  of  the  wave.  In  addition,  there  is 
present  a  region  of  fluid  recirculation  under  the  wave 
when  viewed  in  a  coordinate  system  moving  with  the 


wave.  Thus,  in  the  presence  of  large  waves  the  convec¬ 
tion  and  circulation  are  expected  to  significantly  en¬ 
hance  the  mass  transfer  rale. 

At  Re  -  230,  Fig.  8  illustrates  the  effect  of  the  wave 
motion  even  when  the  waves  are  not  large.  In  general, 
peaks  in  the  values  of  <p  appear  near  the  front  and 
back  of  the  wave.  When  the  wave  is  asymmetric  there 
is  a  sharp  decrease  in  4>  below  the  peak  of  the  wave, 
indicative  of  a  low  concentration  in  the  recirculation 
region.  By  contrast,  the  large  waves  associated  with 
the  higher  flow  rate  in  F;g.  9  all  show  large  peaks  in 
the  solute  distribution  coefficient  in  the  front  region  of 
the  wave.  The  amount  of  mass  captured  by  these  large 
waves  (proportional  to  <f> )  is  2-5  times  larger  than  that 


Normalized  Film  Thickness,  Solute  Distribution  Coefficient 

Fig.  10.  Probability  density  functions  of  normalized  film  thickness  and  solute  distribution  coefficient  for 
Re  —  174  and  a  development  length  of  1.3  m. 
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Fig.  II.  Probability  density  functions  of  normalized  film  thickness  and  solute  distribution  coefficient  for 
Re  =  805  and  a  development  length  of  2  4  m. 
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for  the  small-wave  system.  These  data  seem  to  sup¬ 
port  the  qualitative  argument  that  the  large  waves 
control  the  mass  transfer  process  because  of  the  large 
secondary  convections  induced  by  these  waves. 

Standard  techniques  were  used  to  compute  the 
probability  distributions  of  h  and  (j>  (Bendat  and 
Piersol,  1971)  and  their  spectra  (Priestley,  1987). 
Figures  10  and  1 1  compare  the  pdfs  for  high  and 
low  Reynolds  numbers.  The  large  tail  in  the  distribu¬ 
tion  of  film  thickness  for  Re  =  805  reflects  the  pre¬ 
sence  of  large  waves.  The  shift  in  the  probability 
distribution  of  4>  to  the  region  of  the  larger  waves 
shows  again  that  these  waves  of  larger  amplitude 
control  the  mass  transfer  process. 


Figures  1 2  and  1 3  show  the  auto  spectra  of  h  and  <p 
and  the  coherence  function  between  the  two  for  the 
two  Reynolds  numbers.  The  coherence  function  in  the 
frequency  domain  is  the  analog  of  the  cross-correla¬ 
tion  function  in  the  time  domain.  For  these  com¬ 
putations  the  Bartlett-Priestley  spectral  window  was 
used  with  a  frequency  resolution  0.48  Hz.  The  sample 
size  was  such  that  the  standard  random  error  of  12% 
of  the  spectrum  can  be  expected  at  any  frequency. 

At  Re  -  174  (Fig.  12),  the  spectrum  of  the  film 
thickness  peaks  at  about  8  Hz;  however,  considerable 
energy  is  displayed  due  to  the  smaller  waves  at  higher 
frequencies.  The  spectrum  of  F  is  quite  flat,  suggesting 
contributions  to  transfer  from  the  full  range  of  wave 


Fig.  12.  Power  spectral  density  and  coherence  functions  of  normalized  film  thickness  and  solute  distribu¬ 
tion  coefficient  for  Re  =  174  and  a  development  length  of  1.3  m. 


Fig.  13.  Power  spectral  density  and  coherence  functions  of  normalized  film  thickness  and  solute  distribu¬ 
tion  coefficient  for  Re  =  805  and  a  development  length  of  2.4  m. 
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frequencies.  However,  Fig.  13  shows  that  at  Re  =  805 
the  dominant  frequencies  for  both  the  wave  motion 
and  the  mass  transfer  overlap  in  the  range  of  the  large- 
wave  frequencies  with  a  much  smaller  amount  of 
spectral  energy  in  the  region  associated  with  the  small 
waves.  Again,  this  suggests  that  at  higher  flow  rates 
the  large  waves  dominate  the  process  of  mass  transfer. 

SUMMARY  AND  CONCLUSIONS 
Falling  liquid  films  are  covered  with  a  random 
array  of  small  and  large  waves  interacting  in  a  com¬ 
plex  fashion.  The  experimental  technique  developed 
for  this  study  measured  the  instantaneous  local  solute 
concentration  averaged  over  the  thickness  of  the 
liquid  film,  simultaneously  with  a  measurement  of  the 
local  film  thickness.  Data  were  collected  over  a  range 
of  Reynolds  numbers  and  development  lengths.  An 
analysis  of  these  data  revealed  that  the  large  waves 
control  the  rate  of  transfer  to  the  film.  This  result  is 
consistent  with  recent  theoretical  studies,  which 
indicate  convective  circulation  paths  which  exist 
under  the  large  waves. 
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NOTATION 

A  optical  absorption  of  contaminated  solution 

C  local  contaminant  concentration,  kg  m'3 

<C>  contaminant  concentration  averaged  from 

the  wall  to  the  interface,  kgm'3 
D  molecular  diffusion  coefficient,  m2  s' 1 

g  acceleration  due  to  gravity,  ms'2 

h  local  film  thickness,  mm 

hN  Nusselt  film  thickness  (=  3/4  Re  v1/g)l>1 
I  intensity  of  radiation,  W  m~2 

Pe  Peclet  number  ( =  Re  Sc) 

Q  liquid-film  flow  rate  per  unit  perimeter, 
m2  s'1 

Re  film  Reynolds  number  ( =  4 Q/v) 

Sc  Schmidt  number  ( =  v/D) 

Tx  transmittance  of  solution  at  a  wavelength  k 
u  local  streamwise  velocity,  m" 1 

v  local  velocity  normal  to  the  boundary,  m  ~ 1 

Vw  wave  velocity,  m  s' 1 

x  axial  coordinate  in  lab  frame,  m 

y  coordinate  normal  to  the  boundary,  m 

Greek  letters 

t  molar  absorptivity  of  contaminant,  m4  kg ' 1 

H  liquid  absolute  viscosity,  kg  m  ' 1  s  ' 1 


v  liquid  kinematic  viscosity,  m2  s' 1 

p  liquid  density,  kg  m  ' 3 

r*.  wall  shear  stress,  N  m'2 
4>  solute  distribution  coefficient  as  defined  in 
eq.  (4) 
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Abstract — From  experimental  measurements  of  a  free  falling  liquid  film  at  Re  =  880,  four  representative 
large,  evolving  or  interacting  waves  were  selected  for  computational  domains  in  which  the  Navier-Siokes 
equations  were  numerically  solved.  The  algorithm  computed  velocity  and  pressure  fields  within  each  wave, 
as  well  as  the  shape  necessary  to  match  experimental  wall  shear  stress  data.  Results  show  interaction  effects 
significantly  modify  flow  fields,  compared  to  large  solitary  waves.  Waves  having  two  peaks  had  two  closed 
recirculation  regions,  with  a  mixing  layer  separating  them.  The  size  of  the  recirculation  regions  was 
dependent  on  the  extent  of  separation  of  the  wave  peaks.  As  with  solitary  waves,  strong  streamwise 
accelerations  exist,  with  both  location  and  magnitude  varying  with  shape  and  evolutionary  character  of 
the  wave.  Heat  and  mass  transfer  rates  must  be  enhanced  by  these  flow  properties,  which  a-e  shown  to 
have  a  complicated  dependence  on  wave  structure.  Examination  of  the  flow  fields  suggests  parabolic 
streamwise  velocity  profiles  are  generally  deficient,  explaining  shortcomings  experienced  by  hydrodynamic 
models  based  on  such  simple  velocity  profiles. 

Key  Words:  fal'ing  liquid  films,  numerical  methods 


INTRODUCTION 

Investigations  of  hydrodynamic  and  transport  properties  of  thin  falling  liquid  films  remain  a 
fertile  research  area  today.  Thin  liquid  films  are  encountered  in  a  wide  variety  of  industrial 
process  equipment,  including  wetted  wall  absorbers,  falling  film  chemical  reactors,  condensers  and 
vertical  tube  evaporators.  At  flowrates  of  industrial  interest,  falling  films  (even  in  the  absence  of 
gas  flow)  evolve  to  a  highly  irregular  wavy  interface  which  is  generally  considered  quasi-stationary. 
Figure  1  displays  a  short  time  trace  of  such  a  falling  film,  10,000  mean  film  thicknesses  below  the 
inlet.  The  film  surface  is  covered  by  a  complex  array  of  large  and  small  waves  moving  over  a 
substrate  which  is  less  than  the  mean  film  thickness.  The  large  waves,  ranging  in  amplitude  from 
2  to  5  times  the  substrate  thickness,  carry  a  large  fraction  of  the  total  mass  flowing,  and  are 
speculated  to  control  the  rates  of  scalar  transport  (Dukler  1977).  Before  the  heat  or  mass  transfer 
rates  to  such  films  can  be  modeled  it  is  necessary  to  understand  the  velocity  distributions  which 
exist  within  these  waves,  as  well  as  the  evolution  of  the  interface.  The  present  work  focuses  on  these 
questions. 

Previous  modeling  efforts  generally  were  limited  to  either  single,  non-interacting  (solitary)  waves 
of  various  thicknesses,  or  the  intersection  of  small  waves.  While  large,  solitary  waves  and  small 
ripples  on  the  substrate  symbolize  asymptotic  behavior  of  the  flow,  examination  of  film  thickness 
measurements,  as  in  figure  1,  shows  these  limiting  cases  are  not  representative  of  the  flow. 
Numerical  simulation  of  the  isolated,  large  waves  at  the  flow  conditions  of  figure  1  (Wasden  & 
Dukler  1989)  suggests  transport  through  the  film  is  enhanced  by  the  interaction  of  large  wave  peaks 
with  the  relatively  slowly  moving  substrate.  This  same  study  determined  that  the  seemingly  slow 
evolution  of  these  waves  is  responsible  for  significant  deviations  in  the  flow  field  from  those 
determined  for  solitary  waves.  The  more  complicated  case  of  rapidly  evolving  or  interacting  waves 
integrates  large  wave  interactions  with  the  surrounding  substrate  and  the  potentially  more 
important  effects  of  wave  interactions. 

Evolution  of  the  large  waves  far  from  the  inlet  may  be  regarded  as  processes  of  coalescence  or 
splitting,  and  result  from  some  type  of  flow  instability.  Large  waves  do  not  grow  without  limit, 
but  split  into  daughter  waves  when  sufficiently  perturbed.  Many  large  waves  overrun  smaller, 
slower  moving  waves,  sometimes  incorporating  the  smaller  wave,  while  often  passing  over  them 
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without  significant  mass  addition.  The  wide  array  of  large  wave  behavior  illustrates  the  need  for 
local  measurements  to  determine  relative  effects  of  wave  evolution  and  large  wave  induced 
convection  on  transport  enhancement. 

Reliable  experimental  measurements  of  the  velocity  distribution  in  the  films  is  exceedingly 
difficult  due  to  the  extremely  small  film  thickness  («  1  mm),  very  short  passage  time  of  each  wave 
(*60  ms)  and  the  random  location  of  the  interface,  as  shown  in  figure  1.  Non-intrusive  methods, 
such  as  LDA,  do  not  provide  sufficiently  fine  resolution  to  investigate  velocity  profiles.  Thus, 
experimental  measurements  of  hydrodynamic  variables  appear  limited  to  the  time  variation  of  wall 
shear  stress  and  film  thickness.  Correspondingly,  analytical  models  have  been  developed  in  the 
absence  of  hard  data  on  the  true  flow  conditions  which  appear  to  exist  in  the  waves. 

Most  analytical  models  of  both  single  and  interacting  waves  extend  the  concepts  advanced  by 
Kapitza  (1964),  based  on  the  use  of  a  parabolic  velocity  profile  and  assuming  that  the  streamwise 
hydrodynamic  variables  scale  with  the  wavelength.  In  1972,  in  examining  various  models  developed 
to  that  date,  Dukler  concluded  that  all  failed  to  accurately  represent  any  measured  characteristics 
of  the  wave  except  at  Reynolds  numbers  well  below  those  of  industrial  interest. 

Modeling  the  wavy  film  flow  by  a  direct  solution  of  the  Navier-Stokes  equations  is  hampered 
by  numerical  stiffness  imposed  by  the  stress-free  interface;  as  a  result,  convergence  is  difficult  except 
at  the  lowest  flowrates.  Previous  numerical  modeling  has  focused  solely  on  non-interacting  large 
waves.  Bach  &  Villadsen  (1984)  explored  the  application  of  a  finite  element  scheme  to  the  unsteady 
problem  of  solitary  waves  developing  from  initial  perturbations  on  the  smooth  film  for  Reynolds 
numbers  up  to  100.  The  film  Reynolds  number  is  defined  as  Re  =  4Q/v,  where  Q  is  the  mass 
flowrate  per  unit  perimeter  and  v  is  the  kinematic  fluid  viscosity.  Their  work  predicted  that  the  flow 
far  from  the  inlet  would  consist  of  solitary  waves  having  one  general  shape,  a  condition  contrary 
to  experimental  fact  even  at  film  Reynolds  numbers  as  low  as  1.  Kheshgi  &  Scriven  (1987)  applied 
a  finite  element  technique  to  a  problem  with  periodic  boundary  conditions  in  the  flow  direction, 
and  verified  the  evolution  of  infinitesimal  disturbances  as  predicted  by  Orr-Sommerfeld  analyses. 
Their  work  was  limited  to  low  flowrates,  and  failed  to  generate  waveforms  comparable  to  those 
observed  experimentally  for  fully-developed  flow. 

Recent  simulations  of  solitary  waves  at  a  high  Reynolds  number  (880)  proceeded  by  solvmg  the 
Navier-Stokes  equations  in  a  partially  determined  flow  domain  (Wasden  &  Dukler  1989).  This 
work  showed  that  isolated  wave  velocity  is  strongly  dependent  on  wave  shape,  and  provided 
evidence  that  numerous  streamwise  length  scales  existed  in  the  flow.  Further,  it  was  determined 
that  effects  of  wave  evolution  are  important  near  the  solid  boundary,  shifting  the  maximum  wall 
shear  stress  in  front  of  the  film  thickness  peak.  The  use  of  a  parabolic  streamwise  velocity  profile 
to  describe  the  flow  in  large  waves  was  shown  to  be  inappropriate  over  a  large  portion  of  these 
waves,  suggesting  analytical  models  based  on  such  velocity  profiles  are  fundamentally  inadequate. 

At  present,  no  suitable  models  for  evolving  waves  exist.  In  the  absence  of  such  models,  and 
experimental  methods  for  measuring  velocity  profiles  in  thin  wavy  films,  a  series  of  numerical 
experiments  was  performed.  The  experiments  propose  to  illuminate  the  subject  of  transport 
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Figure  1 .  Film  thickness  time  trace:  Re  =  880. 
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enhancement  due  to  waves  by  solving  for  the  hydrodynamics  within  typical  experimentally 
determined  waveforms.  For  a  film  Reynolds  number  of  880.  wave  shapes  and  wall  shear  stress  data 
were  measured  in  our  laboratory.  Four  representative  large,  interacting  o;  evolving  waveforms  were 
converted  for  use  as  domains  for  a  finite-difference  code  developed  specifically  for  free  surface 
problems.  The  use  of  experimentally  measured  waveforms  ;s  a  novel  concept,  providing  a  somewhat 
simpler  numerical  task  while  insuring  the  results  will  net  represent  isolated  or  idealized  cases  of 
film  flow.  The  results  of  these  computations  demonstrate  the  transpor*  enhancement  properties 
peculiar  to  evolving  large  waves,  and  present  data  useful  for  future  model  development. 


EXPERIMENTAL  PROCEDURE 

Flow  loop 

For  fully-developed  wavy  film  flow,  film  thickness  and  wai!  shear  stress  data  were  collected  in 
a  50.8  mm  i.d.  vertical  test  section  in  a  flow  loop  described  by  Zabaras  et  al.  (1986).  After  being 
pumped  through  a  calibrated  rotameter,  the  fluid  entered  the  column  through  an  annulus  whose 
inner  wall  consists  of  a  stainless  steel  porous  sinter,  having  100  pm  pore  size.  Combined  with  careful 
leveling  of  the  column  prior  to  data  collection,  this  entry  section  ensured  minimal  deviations  from 
axisymmetric  flow,  and  produced  a  smooth  inlet  flow.  The  measuring  station  was  located  3. 1  m 
(roughly  10,000  mean  film  thicknesses)  below  this  entry  section. 

Measuring  station  and  measurement  techniques 

The  measuring  station,  shown  in  figure  2,  is  patterned  after  that  described  by  Zabaras  et  al. 
(1986).  The  removable  section  allows  simultaneous  measurement  of  film  thickness  and  wall  shear 
tress  at  one  location,  as  well  as  film  thickness  at  another.  The  station  was  constructed  of  the  same 
Nterial  as  the  flow  loop,  and  was  carefully  machined  to  ensure  a  smooth  transition  to  the  station. 
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Film  thickness  probes  consisted  of  twin  parallel  platinum  13%  rhodium  wires.  O.OSmmdia 
spaced  2.5  mm  apa.^,  which  penetrated  the  flow.  Described  in  detail  by  Brown  a  al.  (1978),  a  linear 
relation  exists  between  the  resistance  of  the  film  between  the  wires  and  the  film  thickness. 
Calibration  proceeded  by  setting  the  measuring  station  horizontal,  blocking  the  ends  and 
introducing  different  fluid  levels,  determined  to  within  lO^tm  using  a  cathetometer,  followed  by 
measurement  of  the  resulting  resistance.  Downstream  electronics  for  convert  ng  this  resistance  to 
a  d.c.  voltage  signal  are  described  elsewhere  (Zabaras  el  al.  1986).  Conductance  of  the  fluid  was 
monitored  closely  at  all  times  during  the  calibration  and  data  collection  procedures  to  insure  proper 
correction  of  any  thermally  induced  conductance  drift. 

Wall  shear  stress  measurements  were  based  on  the  electrochemical  mar ,  transfer  method 
described  by  Hanratty  &  Campbell  (1983).  For  the  present  series  of  measurements,  the 
iodine/tri-iodide  system  was  chosen.  The  working  solution  contained  0. 1  M  KI  and  0.004  M  I2(s) 
in  demineralized  water,  and  was  replaced  every  2  h  to  minimize  errors  due  to  iodine  evaporation. 
A  dry  nitrogen  atmosphere  was  used  in  the  flow  loop  to  minirr'ze  oxygen  saturation  of  the 
solution.  Fluid  properties  at  25  °C  are:  density,  1010  kg/m’;  absolute  viscosity,  8.5  x  10~4kg/ms; 
and  surface  tension,  7.12  x  !0~:N/m.  The  cathode  for  this  system  consists  of  a  flush-mounted 
strip  of  platinum  foil,  0.075  mm  (in  the  flow  direction)  x  1  mm  wide,  embedded  in  Plexiglas 
to  insure  electrical  isolation.  By  measuring  the  current  produced  by  an  electrochemical  reaction  at 
the  surface  of  the  cathode,  the  wall  shear  stress  at  that  location  is  determined.  For  the  redox 
reaction 


If  +2e  -*  31“  (cathode), 

31  ~  -♦  If  +  2e  (anode), 

a  concentration  boundary  layer  develops  on  the  cathode  surface,  which  is  polarized  at  —0.8  VDC 
to  insure  the  concentration  approaches  zero.  For  the  iodine  system,  the  range  of  polarization 
voltage  is  quite  broad,  insuring  large  increases  in  flowrate  will  not  deplete  the  electron  source  at 
the  cathode.  Details  concerning  the  downstream  electronics  and  calibration  associated  with  this 
measurement  technique  are  to  be  found  elsewhere  (Zabaras  el  al.  1986). 

It  is  now  recognized  (Mao  &  Hanratty  1985)  that  the  response  of  the  •  lectrochemical  probe  is 
highly  dependent  of  the  nature  of  the  “input”  wall  shear  stress.  For  the  ionic  system  employed  in 
this  study,  errors  in  both  phase  and  magnitude  are  expected  to  be  small  due  to  the  large  (103  s'1) 
mean  velocity  gradient,  small  cathode  surface  area  and  large  Schmidt  number  (v/D)  of  the  fluid 
(^780).  The  relationship  given  by  Hanratty  &  Campbell  (1 983''  between  cathode  current  and  wall 
shear  stress  was  used  in  this  study,  as  the  frequencies  in  the  date  were  sufficiently  low  to  allow  the 
use  of  a  quasi-steady  analysis. 


Data  collection,  processing  and  analysis 

Voltage  signals  from  two  film  thickness  probes  and  the  wall  shear  stress  probe  were  first 
low-pass  filtered  at  1  kHz,  then  fed  to  a  microcomputer-based  A/D  converter.  Each  signal  was 
digitized  at  1  kHz  by  a  data  translation  12-bit  A/D  converter  installed  in  a  DEC  Micro  11/73 
microcomputer.  The  data  set  comprised  1  min  of  data,  and  was  stored  on  the  system  Winchester 
disc  prior  to  applying  calibration  curves  and  writing  the  data  to  magnetic  tape  for  further  analysis. 
Digitization  and  collection  errors  are  expected  to  be  negligible  for  all  d:.  j,  while  calibration  errors 
for  the  film  thickness  measurement  are  expected  to  be  <3%.  Errors  innerent  in  applying 
steady-state  wall  shear  stress  calibration  curves  depend  on  the  nature  of  the  input  signal,  requiring 
separate  examination  of  individual  results.  Zabaras  (1985)  reports  estimated  errors  of  <7%  for 
this  technique. 

Film  thickness  and  wall  shear  stress  data  was  examined  to  locate  typical  interacting  or  evolving 
waveforms.  Due  to  the  limited  amount  of  data  provided  by  the  two  film  thickness  probes, 
determining  the  evolutionary  character  of  waveforms  was  difficult.  However,  several  types  of  wave 
shapes  appeared  frequently,  although  their  individual  dimensions  varied  considerably.  Four 
waveforms  were  chosen  as  representative  of  the  large  wave  structures,  each  having  peak  thicknesses 
greater  than  twice  the  surrounding  substrate  thickness.  In  addition,  each  was  preceded  by  a 
reasonably  smooth  substrate,  shown  by  wall  shear  stress  measurements  to  be  free  of  acceleration. 
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Each  measured  film  thickness  profile  was  converted  from  the  time  domain  into  a  spatial  domain 
for  use  with  the  numerical  algorithm.  Wave  peak  passage  times  between  probes  were  used  to 
estimate  “wave  velocities" — these  values  provided  initial  estimates  for  the  numerical  computations. 


NUMERICAL  METHOD 

Solution  of  free  boundary  problems  requires  methods  for  both  hydrodynamic  calculations  and 
shape  determination.  The  velocity  and  pressure  fields  within  the  wave  were  determined  by  solving 
the  Navier-Stokes  equations  in  primitive  variable  form.  For  a  film  Reynolds  number  of  880,  the 
wave  thickness  generally  was  <  1%  of  the  pipe  radius  and  therefore,  a  two-dimensional  cartesian 
coordinate  system  was  chosen.  The  transformation  of  time  traces  of  film  thickness  to  this 
coordinate  system  comprised  the  shape  determination  portion  of  the  overall  algorithm.  The 
common  method  of  computing  the  free  surface  position  is  to  compute  film  thickness  values  at  fixed 
streamwise  locations.  The  present  work  inverts  this  procedure;  for  a  measured  sequence  of  film 
thickness  values,  streamwise  locations  associated  with  each  value  are  determined  such  that  the 
resulting  shape  and  flow  field  within  satisfies  all  interfacial  boundary  conditions. 

To  develop  the  methodology  for  treating  evolving  films,  waves  were  initially  modeled  as  though 
their  shape  remained  constant  with  time — these  waves  are  termed  solitary.  A  new  streamwise 
coordinate,  z,  is  fixed  on  the  front  of  the  wave  and  extends  in  the  opposite  direction  of  gravity. 
The  film  thickness  profile  in  the  time  domain,  h(t,),  was  converted  to  the  length  domain,  h(z,), 
through  the  transformation 

*i**o+  [1} 

for  /  ranging  from  1  to  the  number  of  film  thickness  points  in  the  isolated  wave.  Through  this 
transformation,  the  wave  profile  was  “stretched”  for  use  as  a  computational  domain,  and  time  was 
removed  from  the  problem.  In  this  coordinate  system,  the  wave  remains  fixed,  and  the  wall  moves 
upward  at  a  constant  speed  given  by  Fw,  the  wave  velocity  for  the  solitary  wave. 

It  is  useful  to  define  a  new  streamwise  velocity  component, 


«(z,y)=  -u'(x,y,t)  +  Vm, 


[2] 


where  u\x,  y,  t)  is  the  streamwise  velocity  in  a  coordinate  system  fixed  on  the  wall.  The  governing 
equations  for  this  viscous,  incompressible  and  isothermal  flow  relative  to  the  moving  wave  become 
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where  v  is  the  velocity  in  the  normal  (y)  direction,  P  is  the  pressure,  g  represents  gravitational 
acceleration  and  v  and  p  are  the  kinematic  viscosity  and  density  of  the  fluid,  respectively.  At  the 
stress-free  interface,  y  =  A(z),  tangential  and  normal  stress  balances  require 
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where  a  is  the  surface  tension  coefficient  and  y.  is  the  absolute  viscosity.  At  the  wall,  y  =  0, 

u  =  V„,  t>=0,  [8] 

represent  the  standard  no-slip  and  no-flux  conditions. 

Velocities  at  the  interface  are  related  through  the  kinematic  condition  in  a  moving  frame. 


The  inlet  velocity  profile  is  parabolic,  representing  an  acceleration-free  falling  film,  while  a  sufficient 
and  physically  consistent  outlet  condition  for  a  solitary  wave  requires  a  zero  streamwise  derivative 
for  all  variables.  The  variable  Vm  replaces  h(z)  as  the  final  variable  to  be  iteratively  determined 
in  the  free  surface  problem. 

For  each  wave  profile,  a  unique,  non-uniform  finite  difference  grid  mesh  was  constructed.  The 
mesh  for  a  typical  domain  is  shown  in  figure  3.  The  particular  wave  shape  determined  the  grid 
spacing  used.  Mesh  refinement  continued  until  no  further  change  in  either  the  computed  wave 
velocity  or  wall  shear  stress  profile  was  observed.  Of  particular  importance  was  the  concentration 
of  cells  near  the  front  of  the  wave,  since  the  velocity  fields  change  drastically  in  this  region  due 
to  the  large  interfacial  slope.  In  addition,  large  curvatures  exist  at  each  peak  and  trough.  To  insure 
adequate  resolution  of  capillary  pressure  induced  effects,  grids  were  concentrated  near  the  free 
interface  in  these  regions.  For  most  waves,  1800  cells  of  dimension  SxSy  were  sufficient,  and 
produced  grid  Reynolds  numbers  (Re0i  =  u(z,  y)dx/v,  ReC/  =  jj(z,y)6>>/v]  of  order  1  in  the  y 
direction,  and  ranging  from  l  to  100  in  the  streamwise  direction. 

The  curved  interface  was  accomodated  by  allowing  boundar  ceils  to  be  cut  by  the  boundary, 
h(z),  thus  reducing  their  volume.  This  situation  is  illustrated  in  figure  4.  This  technique  produced 
areas  adjoining  two  boundary  cells,  the  centers  of  which  were  outside  the  computational  domain. 
As  the  stress-free  interface  requires  a  zero  normal  derivative  of  the  velocity  vector  with  respect  to 
the  boundary,  h(z),  these  regions  were  treated  as  inviscid  channels  through  which  all  fluid  leaving 
one  boundary  cell  on  its  shared  side  passed  into  the  neighboring  cell  through  its  respective  shared 
side.  The  total  area  of  these  regions  represents  <0.1%  of  the  total  domain,  and  had  little  effect 
on  the  results. 

The  equation  set  was  solved  on  a  finite  difference  grid  using  a  variant  of  the  TEACH-T  code 
(Gosman  et  al.  1969),  incorporating  the  SIMPLER  pressure/continuity  solution  procedure;  the 
principles  of  this  method  are  described  in  detail  elsewhere  (Patankar  1980).  The  domain  includes 
regions  of  significant  streamwise  variation  in  all  variables,  thus  necessitating  an  accurate  method 
of  discretizing  convective  momentum  terms.  The  simplest  method  of  convective  discretization, 
upwind  differencing,  ensures  a  reasonably  stable  numerical  solution,  but  introduces  numerical 
diffusion  in  regions  of  the  flow  where  streamlines  are  oblique  with  respect  to  the  grid  lines 
(Raithby  1976).  More  importantly,  the  upwind  scheme  lacks  sensitivity  to  cross-stream  diffusion 
and  source  terms  (Leonard  1979),  which  are  of  tremendous  importance  in  the  case  of  a  thin  film. 


Tim*,  ms 

Figure  3.  Sample  finite  difference  grid. 
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This  lack  of  sensitivity  diminishes  the  effects  of  the  y  direction  diffusion  as  well  as  the  "t  "  velocity 
in  the  solution  of  the  streamwise  velocity.  These  deficiencies  in  the  upwind  and  hybrid  methods 
require  the  use  of  a  QUICK-based  scheme,  which  improves  acccuracy  by  expanding  the 
number  of  neighboring  points  included  in  interpolated  values  of  velocity.  Based  on  Leonard’s 
(1979)  third-order  accurate  discretization  scheme  QUICK,  Pollard  &  Siu  (1982)  developed  the 
QUICK-ER  (Extended  and  Revised)  method  of  discretizing  convective  terms.  The  QUICK-ER 
method  overcomes  stability  problems  inherent  in  the  QUICK  procedure  at  the  expense  of 
slower  convergence,  and  is  considered  the  most  satisfactory  method  of  handling  convective 
momentum  terms  (Huang  el  al.  1985).  For  application  to  non-uniform  grids,  a  new  version  of 
QUICK-ER  was  developed.  This  method  follows  the  spirit  of  the  QUICK-ER  formulation,  but 
includes  locally  variable  weighting  factors  to  account  for  the  non-uniformity  of  the  grid  in  both 
directions.  Although  QUICK-ER  schemes  requires  more  computational  effort  per  iteration  than 
upwinding,  particularly  for  non-uniform  grids,  improvements  in  accuracy  allow  the  use  of  a 
slightly  coarser  grid,  so  total  computational  time  exceeds  that  required  by  the  upwind  method  by 
only  20%. 

The  solution  procedure  began  with  choosing  a  value  for  Vw  and  creating  the  transformed 
domain,  given  by  [1].  The  u  velocity  field  was  set  to  a  parabolic  profile  everywhere,  and  the  v 
velocity  field  was  set  to  zero.  The  pressure  at  each  z  location  was  set  to  the  surface  pressure  due 
to  curvature.  Updated  velocity  and  pressure  fields  within  the  wave  were  then  computed  using 
[3]— [5].  Through  interpolation  for  the  velocity  gradients  in  the  interfacial  shear  stress  balance,  [6], 
streamwise  and  normal  velocities  in  the  interior  of  the  flow  field  were  used  to  derive  an  expres¬ 
sion  for  the  streamwise  surface  velocity.  Coupled  with  the  kinematic  condition  (9],  the  velocities 
on  the  surface  were  known  for  each  iteration.  The  surface  pressure  computed  from  (7]  was  used 
to  determine  the  first  pressure  value  in  the  interior  of  the  domain  through  the  use  of  parabolic 
interpolation  using  the  surface  pressure  and  two  interior  pressures.  With  the  newly  computed 
surface  variables,  the  velocity  and  pressure  fields  were  updated  until  the  sum  of  residuals 
of  mass  and  momentum  (normalized  by  the  inlet  quantities)  over  the  domain  was  <  I0'3. 
In  addition,  the  interfacial  shear  and  normal  stress  balances  were  required  to  be  within  10' 2  Pa 
of  zero. 

The  free  boundary  shape  was  determined  by  examining  the  converged  average  outlet  pressure 
in  the  flat  film  section.  If  the  average  pressure  did  not  approach  zero,  as  required  for  a 
non-accelerating  film  trailing  a  solitary  wave,  a  new  value  of  Fw  was  chosen  and  the  process 
repeated.  The  adjustment  procedure  for  Vw  was  simple — if  the  pressure  in  the  outlet  section  was 
higher  than  zero,  the  wave  (wall)  velocity  was  too  high,  since  the  wall  was  pushing  excess  fluid 
through  the  wave,  and  a  positive  pressure  at  the  outlet  was  opposing  this  extra  fluid  in  an  attempt 
to  satisfy  the  mass  balance  for  the  wave. 

The  numerical  study  of  evolving  or  interacting  waves  retains  much  of  the  methodology  developed 
for  solitary  waves,  differing  only  in  the  wave  shape  determination.  Examination  of  the  experimental 
measurements  of  film  thickness  reveals  the  large  waves  change  rapidly  between  the  upper  and  lower 
film  thickness  probe.  The  waves  appear  to  change  from  solitary  type  waves  by  being  perturbed  by 
locally  variable  mass  and  momentum  sources,  physically  recognized  as  small  waves.  Incorporating 
this  unsteady  effect  is  accomplished  through  the  use  of  a  locally  constant  stretching  parameter,  as 
opposed  to  the  globally  constant  value  used  for  the  classical  solitary  wave.  The  domain 
transformation  for  this  case  is  given  by 

z,  =  z0+  l10l 

where  z  is  the  streamwise  coordinate  fixed  on  the  wave  front.  In  general,  this  pseudo  wave 
velocity  is 

Kw.,-Kw[l  +e(z,)],  (11] 

where  e(z ,)  is  an  iteratively  determined  local  stretching  variable  and  Vw  represents  the  wave  velocity 
associated  with  the  substrate.  The  solitary  wave  case  is  recovered  by  setting  e(z,)  =  0  Vi.  As  before, 
we  define  a  new  streamwise  velocity  component  as 


u(z,y)=  —u'(x,y,  t)  +  K„(z,). 


[12] 
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and  continuity  equations  [3]  and  [5],  which  become 
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The  normal  stress  balance,  [16],  was  found  to  be  insensitive  to  the  additional  term,  dKw/dz,  and 
the  term  was  removed  from  the  formulation.  The  solution  procedure  for  the  pseudo-unsteady 
equations  [13],  [4]  and  [14],  subject  to  [15]  and  [7]-[9],  is  identical  to  that  of  the  solitary  wave  with 
the  exception  that  now  a  profile  of  Vw ,  must  be  specified  instead  of  a  single  value.  When  the  velocity 
and  pressure  fields  have  converged  for  a  given  set  of  VmJ,  the  wave  shape  is  adjusted  through  e(z,) 
to  meet  two  criteria.  The  baseline  wave  velocity,  Vv,  is  adjusted  such  that  the  average  pressure  in 
the  flat  outlet  section  approaches  zero,  as  before.  The  computed  wall  shear  stress  profile  is  then 
compared  to  the  experimental  profile,  and  adjustments  made  to  e(z,)  to  correct  deviations.  Upon 
arriving  at  the  correct  distribution  of  e(z,),  the  solution  includes  the  velocity  and  pressure  fields 
within  the  wave,  as  well  as  the  relative  extents  of  its  evolution  throughout  the  wave.  Each  solution 
is  unique,  and  the  wall  shear  stress  matching  procedure  insures  the  accuracy  of  the  velocity  fields. 

The  procedure  developed  for  the  solitary  waves  required  an  average  of  300  iterations  of  the 
velocity  and  pressure  fields  to  converge,  with  an  under-relaxation  factor  of  one-half  used  for  all 
variables.  Between  4  and  8  adjustments  to  the  solitary  wave  velocity  were  required  to  produce  a 
flow  with  an  average  outlet  pressure  <  10-J  Pa.  For  the  quasi-unsteady  case,  roughly  500  iterations 
were  required  to  achieve  convergence  of  the  velocity  and  pressure  fields,  due  to  the  extra  stiffness 
imposed  by  the  multiple  peaks  and  valleys.  Adjustment  of  the  variable  wave  velocity  to  match  the 
wall  shear  stress  data  took  anywhere  from  20  to  40  iterations. 

The  program  was  coded  in  FORTRAN  77,  and  required  2  mbyte  of  task  space.  Execution  times 
for  convergence  of  the  velocity  and  pressure  fields  for  a  given  wave  shape  were  approx.  5  CPU  h 
on  a  VAX  11-750. 
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This  study  concentrated  on  waves  developing  from  or  interacting  with  isolated  waves,  identified 
by  comparing  film  thickness  measurements  from  both  upper  and  lower  probes.  Examination  of  the 
experimental  film  thickness  data  shows  that  many  large  waves  fall  into  three  categories.  While 
nearly  solitary  waves  ex;st,  few  are  free  from  small  wave  induced  ripples,  an  example  of  which  is 
shown  in  figure  5.  As  the  wave  travels  from  the  upper  to  the  lower  probe,  the  small  hump  on  the 
wave  tail  grows,  while  the  wave  front  retains  its  structure;  we  classify  this  wave  as  type  A.  Large 
waves  appear  susceptible  to  splitting  when  perturbed,  as  illustrated  in  figure  6.  These  waves, 
classified  as  type  B,  cover  a  large  portion  of  the  interface,  and  suggest  that  modeling  falling  films 
far  from  the  inlet  as  steady  processes  neglects  important  dynamics.  The  interaction  of  two  large 
waves  of  unequal  size  is  shown  in  the  wave  trace  in  figure  7.  Although  significantly  smaller,  the 
trailing  wave  seen  in  the  upper  trace  either  passes  through  the  larger  wave  or  accumulates  a  larger 
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Figure  5.  Evolution  of  a  nearly  solitary  wave:  type  A.  Figure  6.  Evolution  of  a  large  wave  via  splitting:  type  B. 


portion  of  its  mass,  leading  to  a  situation  classified  as  type  C.  While  each  of  these  cases  characterizes 
a  significant  fraction  of  the  interface,  combinations  of  these  effects  are  also  seen.  Figure  8  illustrates 
the  splitting  of  a  large  wave  coupled  with  the  interaction  of  the  large  wave  with  a  small  wave 
preceding  it.  This  type  D  wave  is  not  representative  of  all  complicated  waves,  but  serves  as  an 
example  of  the  combination  of  evolution  and  interaction  effects. 

Description  of  large  waves  is  generally  limited  to  peak  and  substrate  thicknesses  and  wave¬ 
lengths.  Experimental  data  is  limited  to  two  time  traces  of  film  thickness,  which  makes  quantitative 
evaluation  of  these  variables  difficult.  It  must  suffice  to  characterize  the  extent  of  evolution  by  the 
time  separation  of  individual  peaks,  while  the  relative  mass  of  each  peak  is  estimated  on  the  basis 
of  the  peak  thickness.  In  conjunction  with  the  numerical  predictions  of  flow  fields,  these  data  offer 
new  insight  into  the  hydrodynamic  processes  occurring  in  the  large  waves. 

Streamlines  computed  for  the  type  A  wave  are  shown  in  figure  9.  The  presence  of  a  large 
recirculation  region  is  seen  within  the  peak  of  the  wave  when  viewed  in  a  coordinate  system  fixed 
on  the  wave.  The  interaction  of  this  large  mass  of  fluid  with  the  surrounding  substrate  results  in 
regions  of  strong  streamwise  acceleration  ahead  and  behind  the  region.  These  results  were  also  seen 
in  studies  of  nearly  solitary  waves  with  similar  peak/substrate  thickness  ratios  (Wasden  &  Dukler 
1989).  The  presence  of  the  growing  hump  on  the  wave  tail  increases  the  size  of  the  recirculation 
region,  compared  to  a  solitary  wave,  which  causes  larger  disruptions  in  the  substrate  as  the  wave 
peak  passes  over.  The  wall  shear  stress  profile  for  this  wave  is  shown  in  figure  10,  along  with 
experimental  values.  Agreement  is  excellent  over  most  of  the  wave — the  experimental  values  at  the 
tail  suggest  outside  or  cross-stream  disturbances  are  present.  This  discrepancy  is  not  expected  to 
significantly  alter  the  overall  flow  field,  as  the  region  under  the  peak  is  insensitive  to  downstream 
disturbances  at  Re  =  880.  Matching  the  wall  shear  stress  data  required  stretching  the  wave  shape 
in  transforming  it  from  the  time  domain  to  the  spatial  domain,  as  shown  in  figure  11.  The 
transformation  required  that  the  entire  wave  stretch  in  length,  compared  to  a  purely  solitary 
wave;  the  magnitude  of  e(zf),  [11],  is  shown  in  figure  II.  Of  special  interest  is  the  extra 
stretching  caused  by  the  presence  of  the  hump  on  the  tail — previous  studies  of  nearly  solitary 
waves  (Wasden  &  Dukler  1989)  show  the  evolution  is  confined  to  the  region  very  near  the 
peak,  while  in  this  case,  the  tail  region  is  stretching  as  well.  From  a  substrate  wave  velocity  of 
1. 10 m/s,  the  wave  velocity  Vw/,  [11],  increased  to  1.89  m/s  near  the  peak.  The  average,  or 
Nusselt,  film  velocity  for  this  Re  (Kapitza  1964)  is  0.51  m/s,  showing  the  wave  moves  at  between 
2.1  and  3.7  times  the  average  velocity.  The  combination  of  experimental  observations  and 
numerical  simulations  show  that  the  wave  is  in  the  early  stages  of  splitting,  and  that  this  effect, 
although  seemingly  slight  upon  examining  the  film  thickness  measurements,  significantly  affects 
the  flow  field. 

Flow  within  a  wave  with  the  peak  splitting  into  daughter  waves  is  shown  in  figure  12.  The  flow 
beneath  each  peak  resembles  that  seen  in  the  nearly  solitary  wave  A.  The  recirculating  region 
encloses  the  entire  area  under  the  twin  peaks,  and  again  causes  large  accelerations  near  the  front 
and  rear  of  the  wave.  The  remnants  of  the  original  solitary  wave  recirculation  region  are  seen  at 


Film  Thickness,  mm 


366 


F.  K.  WASDEN  »nd  A  C  DUKLER 


Tims,  ms  Time,  ms 

Figure  7.  Interaction  of  large  waves:  type  C.  Figure  8.  Evolution  and  interaction  of  large  waves:  type  D. 
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m  the  front  of  the  wave,  with  the  newly  formed  rear  peak  having  a  smaller  recirculation  region,  even  prev 

though  it  has  a  roughly  equivalent  peak  thickness.  This  observation  demonstrates  the  danger  of  The 

■  classifying  multiple  waves  based  solely  on  the  easily  measured  parameters  such  as  peak  thick-  The 

|  nesses  and  separation  time.  Comparison  of  computed  and  measured  wall  shear  stress  is  shown  in  Nus 

figure  13  to  be  excellent.  The  peak  in  wall  shear  stress  preceding  the  first  peak  in  film  thickness  T1 

H  is  similar  to  that  observed  in  an  isolated  wave,  showing  how  the  evolving  wave  retains  part  of  its  with 

■  original  traits.  The  second  shear  stress  peak  is  a  further  indication  of  the  emergence  of  a  new  peak.  wavt 

The  domain  and  stretching  parameter  variation  for  this  wave  are  shown  in  figure  14.  The  stretching  com] 

parameter  was  varied  such  that  the  peaks  were  evolving  most  rapidly,  with  the  trough  moving  accel 

■  slower  than  either  peak.  As  expected,  the  second  peak  was  moving  slightly  (» 10%)  slower  than  flow. 

■  the  first.  The  substrate  wave  velocity  was  1.32  m/s,  yielding  Vm  t  from  1.32  to  1.91  m/s  under  the  awaj 

peaks.  Compared  to  the  Nusselt  velocity,  the  wave  moves  roughly  2.6-3.8  times  faster.  inter 

Mj  Flow  in  an  interacting  wave  sequence  is  shown  in  figure  15.  The  two  waves  appear  to  have  of  tb 

||  independent  recirculation  regions,  with  the  first  preceded  by  the  characteristic  acceleration  region.  evoh 

The  size  of  the  recirculation  zones  appears  to  be  related  to  the  thicknesses  of  the  individual  peaks,  altho 

—  in  contrast  to  the  evolving  case,  type  B.  This  may  result  from  the  extra  separation  that  exists  in  data 

■  this  case,  causing  speculation  that  the  type  B  wave  would  evolve  into  a  wave  similar  to  type  C.  is  exc 

™  Shear  stress  profiles  for  the  interacting  case  are  displayed  in  figure  16.  The  agreement  between  effect 

computed  and  measured  shear  stress  is  excellent,  and  shows  two  peaks  in  shear  stress,  indicating  stress 

S  the  wave  is  evolving  into  two  nearly  solitary  waves.  As  the  wave  peaks  separate  further,  the  second  This 

■  wall  shear  stress  peak  is  closer  in  magnitude  to  that  expected  of  an  isolated  wave,  additional  subst 

evidence  of  future  disintegration  of  the  original  wave.  The  transformed  domain  and  stretching  and  t 

■j  parameter  variation  are  shown  in  figure  17.  Both  peaks  were  evolving  at  a  nearly  identical  rate,  no  si. 

■  suggesting  the  effect  of  the  interactions  on  the  natural  evolution  of  the  wave  was  slight.  As  was 
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Figure  9.  Streamline  map:  wave  A. 
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Figure  10.  Shear  stress  comparison:  wave  A. 
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Figure  II.  Experimental  and  computed  wave  shape:  wave  A. 


previously  the  case,  the  trough  between  the  waves  was  stretching  more  slowly  than  the  peak  regions. 
The  substrate  wave  velocity  was  1.12  m/s,  closer  to  the  solitary  wave  (A)  than  the  splitting  wave  (B). 
The  wave  velocity,  Kw  i,  varied  from  1.12  to  1.7  m/s  under  the  peaks,  roughly  2.1-3.3  times  the 
Nusselt  velocity. 

The  most  complicated  waveform  studied  included  the  effects  of  evolution  and  interaction.  Flow 
within  this  wave,  classified  as  type  D,  is  shown  in  figure  18.  Apparently  the  interaction  of  the  large 
wave  structure  with  the  small  wave  near  the  front  has  little  effect  on  the  overall  flow  pattern, 
compared  to  what  occurs  in  simple  evolving  waves  of  this  size  (see  figure  1 2).  A  region  of  moderate 
acceleration  exists  at  the  front  of  the  small  forerunner  wave,  followed  by  a  region  of  nearly  parallel 
flow.  This  flow  profile  suggests  the  effect  of  the  initial  interaction  is  limited  to  accelerating  fluid 
away  from  the  wall,  and  effectively  changing  the  substrate  thickness  with  which  the  large  wave 
interacts.  Further  downstream,  the  large  acceleration  region  is  again  evident,  due  to  the  interaction 
of  the  large  recirculation  region  with  the  slowly  moving  substrate.  As  the  wave  has  only  begun  to 
evolve,  the  recirculation  region  is  reminiscent  of  that  associated  with  a  stretched  isolated  wave, 
although  the  area  near  the  trough  of  the  wave  implies  future  separation  of  the  waves.  Shear  stress 
data  for  this  wave  is  shown  in  figure  19;  agreement  between  the  measured  and  computed  values 
is  excellent  over  the  entire  wave.  The  shear  stress  profile  provides  further  evidence  of  the  dampening 
effect  of  the  interaction  of  the  large  wave  with  the  small  forerunner  wave.  The  maximum  shear 
stress  preceding  the  large  wave  peak  is  smaller  (relative  to  pgh)  than  existed  in  the  previous  cases. 
This  may  be  due  to  the  relative  difference  in  size  between  the  wave  peak  and  the  local 
substrate — previously,  wave  peak/substrate  ratios  were  *  3,  while  the  ratio  between  the  first  peak 
and  the  plateau  preceding  it,  in  this  case,  is  <  2.  Since  the  splitting  of  the  peaks  has  just  begun, 
no  significant  secondary  peak  in  shear  stress  exists.  The  transformed  domain  and  stretching 
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Figure  12.  Streamline  map:  wave  B. 
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Figure  13.  Shear  stress  comparison:  wave  B. 
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Figure  14.  Experimental  and  computed  wave  shape:  wave  B. 


parameter  variation  are  shown  in  figure  20,  and  are  similar  to  cases  B  and  C.  Evolution  of  the 
forerunner  wave  is  significant;  the  wave  is  stretching  out  in  response  to  the  mass  descending  on 
it.  As  before,  the  second  peak  is  moving  slightly  slower  than  the  first  (s=5%);  from  data  shown 
in  case  B,  the  difference  is  expected  to  increase  with  the  extent  of  evolution.  From  a  substrate  wave 
velocity  of  1.24  m/s,  the  wave  velocity  Vw  ,  increased  to  1.89  m/s  under  the  peak,  ranging  from 
2.4  to  3.8  times  the  Nusselt  velocity. 


SUMMARY  AND  CONCLUSIONS 

Interacting  or  evolving  large  waves  comprise  a  large  portion  of  the  interface  of  a  falling  liquid 
film.  Wave  amplitudes  of  2-5  times  the  substrate  thickness  are  common  even  for  moderate 
Reynolds  numbers,  ruling  out  models  based  on  small  perturbation  theory.  From  film  thickness 
measurement  time  traces,  four  representative  large  waves  were  selected  as  comutational  domains 
for  a  numerical  solution  of  the  Navier-Stokes  equations.  Computed  results  included  velocity  and 
pressure  fields,  as  well  as  the  wave  shape  necessary  to  match  shear  stress  measurements.  Values 
of  yw,/  111],  remain  near  those  expected  for  solitary  waves  with  similar  peak  and  substrate 
thicknesses,  as  reported  by  Wasden  &  Dukler  (1989).  The  stretching  technique  developed  for 
this  study  allowed  simple  accomodation  of  local  unsteady  effects  by  decoupling  the  hydrodynamic 
and  shape  determination  problems.  An  interesting  extension  of  the  present  work  would  be  to 
measure  the  film  thickness  at  several  closely  spaced  locations,  and  attempt  to  compute  e(z )  from 
the  data. 

The  bulk  of  the  liquid  in  large  interacting  or  evnHng  waves  is  carried  in  the  region  above  the 
substrate,  and  is  nearly  stationary  in  a  coordinate  system  moving  with  the  wave.  As  the  wave  moves 
rapidly  over  the  slow  substrate,  it  causes  acceleration  of  fluid  from  the  front  into  the  peak.  The 
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Figure  1 5.  Streamline  map:  wave  C. 
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fluid  then  decelerates  as  it  passes  out  of  the  wave  tail,  generating  a  closed  recirculation  region  in 
the  wave  peak.  Waves  with  multiple  peaks  have  multiple  regions  of  recirculation,  with  mixing  zones 
between  them.  These  regions  occur  only  in  multiple-peak  waves,  and  may  be  responsible  for  heat 
and  mass  transfer  enhancement  above  that  due  to  fluid  acceleration  and  circulation.  The  magnitude 
of  the  enhancement  due  solely  to  multiple  wave  interactions  is  unknown,  but  may  be  significant 
compared  to  that  associated  with  isolated  large  waves.  The  computations  were  limited  to  waves 
having  two  peaks,  but  these  results  are  expected  to  extrapolate  well  for  cases  of  more  than  two 


■derate 
cldflbs 
>m0ts 
ty  and 


ov^he 
imj'cs 
ik.  The 


Dtatanca,  mm 

Figure  18.  Streamline  map:  wave  D. 
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Figure  19.  Shear  stress  comparison:  wave  D,  Figure  20.  Experimental  and  computed  wave  shape:  wave  D. 
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peaks.  As  a  majority  of  the  surface  is  covered  by  evolving  or  interacting  waves,  further  efforts 
toward  determining  relative  effects  of  competing  hydrodynamic  processes  is  justified. 

Perhaps  the  most  surprising  result  of  this  study  concerned  the  similarity  between  flow  fields 
resulting  from  the  interaction  of  large  waves  and  those  generated  by  the  splitting  of  a  large  wave 
into  two  waves.  This  similarity  suggests  the  shape  of  the  interface  controls  the  fluid  dynamics 
within,  regardless  of  the  evolutionary  process  from  which  the  wave  resulted.  This  may  simplify 
further  analyses  of  transport  processes,  and  allow  characterization  based  primarily  on  peak  and 
substrate  thicknesses,  as  well  as  separation  times. 

Modeling  large  wave  interactions  will  require  the  use  of  velocity  profiles  capable  of  representing 
the  large  streamwise  accelerations  existing  in  the  flow.  Previous  work  (Wasden  &  Dukler  1989)  has 
shown  a  cubic  polynomial  reasonably  approximates  the  normal  (y)  variation  of  the  streamwise 
velocity  at  all  axial  locations:  its  use  for  modeling  the  evolution  of  large  waves  is  recommended. 
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